††thanks: Dena Monjazebi was with the Department of Science, Ryerson
University, Toronto, ON, Canada e-mail<EMAIL_ADDRESS>
# Estimating Phase Aberration from Noisy Radiofrequency Data of a single frame
of Synthetic Aperture Ultrasound Image

1st Dena Monjazebi  dept. Physics
Ryerson University
Toronto, Canada
<EMAIL_ADDRESS>2nd Yuan Xu  dept. Physics
Ryerson University
Toronto, Canada
<EMAIL_ADDRESS>
###### Abstract

Phase aberration is one of the main contributors to image degradation in
ultrasound imaging. Normalized-Cross-Correlation (NCC) is one of the most
extensively studied techniques to estimate the arrival delay error and the
aberration profile. However, the performance of NCC can be compromised when
the data have a poor signal-to-noise-ratio. Here we propose an iterative
filtered NCC (f-NCC) method to estimate phase aberration in speckle regions
from the noisy radio-frequency data of a single frame of synthetic transmit
aperture (STA) image. First, two-dimensional filters were applied in the
aperture and temporal domains to suppress noise and off-axis signals in the
transmit beamformed data from an arbitrary speckle region. Second, an
iterative method that included the phase aberration in both the transmission
and receiving process was used to estimate the phase aberration. Third, the
estimation were improved by correcting the edge effect with an iterative zero-
padding technique. Last, in the experimental study, we proposed and validated
a formula to estimate the arrival time error caused by a sound speed error.
The results from both the numerical simulations and experimental studies
confirm that the proposed filtered NCC method can estimate the phase
aberration profile correctly from noisy STA RF data in a speckle region. Image
quality in terms of contrast, resolution and image signal-to-noise-ratio was
improved. Potential applications of the proposed method to map local speed of
samples are also discussed.

###### Index Terms:

Phase Aberration, Cross-Correlation, Transmit Beamformed Data, Two-Dimensional
Filter, Aperture and Temporal Domain

##  I Introduction

Phase aberration can compromise the performance of both ultrasound imaging and
ultrasound therapy. Recent studies has shown that phase aberration due to
sound speed heterogeneity is the dominant cause of image degradation resulted
from poor beam-focusing quality[1, 2, 3, 4, 5]. In the presence of speed
heterogeneity, delay calculations solely based on the geometric analysis are
incorrect. These errors will cause poor focusing, reduce contrast and worsen
lateral resolution. With fatty tissue having a different sound-speed compared
to other organs, phase aberration is a major problem in ultrasound imaging of
breast and obese patients.

Phase aberration should be compensated to align the wavefront during
beamforming. Many approaches have been developed to correct the phase
aberration. Some of these methods are: normalized cross-correlation (NCC) [6,
7, 8, 9, 10], maximization of speckle brightness [11], time reversal of
ultrasound field [12] and minimization of the sum of absolute differences
between signals at adjacent elements [13]. Aperture domain model image
reconstruction (ADMIRE) has been proven effective to reduce the acoustic
clutter by decomposing the aperture domain data to differentiate between
clutter and main beam [14].

NCC is a widely studied delay estimation algorithm in ultrasound imaging and
therapy. O’Donnell and Flax have demonstrated that normalized and non-
normalized cross-correlation techniques can be used to obtain accurate
estimations of delay compensation in B-mode imaging [6][10]. Generally, a
segment of radiofrequency (RF) data in a reference received signal was
compared with a segment of RF data in a delayed received signal. A pattern-
matching function was utilized to find the delay that maximized the
resemblance of the two segments [15]. Regarding commonly used algorithms,
normalized and non-normalized NCC gave reasonable results in terms of
precision and computational time [16]. Many papers have used the NCC
aberration profiles information to estimate an average speed of the medium or
a local speed map [17, 18, 19].

Efficacy of NCC studies highly depends on the Signal-to-Noise-Ratio(SNR) of
the data. This requirement introduces a major shortcoming in in vivo studies.
In hard-to-image patients due to the extra intermediate fat layer, organs and
anatomical features are located at a larger depth. With the additional layer
of fat attenuating the signal, RF signals from organs of interest have
inherently lower SNR. Slack et al. have indicated prevalence of obesity in US
as 30.8 %. This Study reported a possible projection of 42 % by 2030 [20]. In
clinical abdominal ultrasound, over-weightness is the biggest limitation to
produce optimal quality images[21]. The requirement for high SNR signals,
introduces a major obstacle to expand the clutter suppression methods into
clinic. With the rising rates of obesity, the need for a successful phase
aberration correction method for low $RF_{SNR}$ in ultrasound visualisation is
increasing.

Various methods with the goal to correct the dilemma of phase aberration
correction method in presence of noise has been developed. Many aperture
domain filters have been investigated to reduce the off-axis signals and noise
[22, 23]. Shin et al. have proposed a method with the focus on phase
aberration correction for RF in presence of of noise [24]. Other methods
include spatial compounding to suppress phase aberration and noise clutters
together[25].

Most phase aberration methods were discussed for B-mode imaging. In this
paper, we focus on synthetic transmit aperture (STA) image due to feasibility
of using the RF data of a single frame of STA image to estimate the phase
aberration iteratively. In STA data acquisition, a semi-spherical wave is sent
by a single transmission element to cover a large imaging area. The echo
signals are recorded by the entire receive aperture. This received data can be
used to reconstruct a low-resolution image. By repeating the transmission
across the entire aperture, a three-dimensional STA data set can be acquired
to reconstruct a high-resolution image. STA imaging can achieve the optimum
transmission and receiving focus across the whole image. In addition, STA data
set has more information than that of B-mode imaging [26].

However, it is known that STA RF data exhibits a smaller $RF_{SNR}$ due to
smaller transmission energy than the conventional B-mode data. In this paper,
we propose a solution to improve the NCC method for noisy RF signals by
filtering the RF data in the aperture-temporal domain. A low-pass aperture-
temporal domain filter was designed to pass the main beam signals and filter
out the acoustic clutter and noise. Previous research indicated off-axis
scattering can be reduced by a aperture-domain filter to improve the phase
aberration estimations in B-mode imaging[22, 23]. Here, f-NCC method was
developed for areas with speckles in STA imaging where $RF_{SNR}$ is
significantly lower than B-mode.

In the method section, we will introduce the low-pass aperture-temporal domain
filter that applies gradual cut-offs in aperture and temporal frequency
domains. The efficacy of this filter to minimize the clutter and off-axis
scattering in the transmit-beamformed (TBF) signals will be investigated. In
the results section, the f-NCC algorithm will be studied in simulation and
experimental level. Alternative potential of f-NCC methodology and future work
will be explained in the discussion section. Finally the results and
conclusions will be summarized in the conclusion section. The appendices
include an introduction to NCC, zero-padding to alleviate the edge effect, and
the derivation of the phase aberration profile in a average speed error
scenario.

##  II METHOD

###  II-A Iterative estimation of phase aberration in STA

In this study, we utilized an NCC-based iterative method to estimate the phase
aberration [13] in Synthetic Transmit Aperture (STA) ultrasound imaging since
it does not require multiple rounds of data acquisition. Let us first define a
complete set of STA RF data: $P_{ij}(t)$, which was the signal at time $t$
from the $i_{th}$ ($i=1:N$) receive element upon the activation of the
$j_{th}$ ($j=1:M$) transmission element. We used $g_{i,j}(x_{f},z_{f})$ to
represent the focusing or geometrical delay corresponding to the focal point
at $x_{f}$ and $z_{f}$. The geometrical delay was defined as the travel time
from transmit element $j$ to the focal point and back to the receive element
$i$. We used $s_{i}(t)$ to represent the transmission beamformed (TBF) RF
signals. This represents the RF signals received by the $i_{th}$ ($i=1:N$)
receive element when the transmitted beam was focused at the focal point
($x_{f}$, $z_{f}$).

TBF signal $s_{i}(t)$ was obtained by adopting the Delay-and-Sum (DAS)
algorithm: summing over all transmit elements of the RF signals from the
$i_{th}$ receive element (Eq. 1) after applying the geometrical delay $g_{ij}$
and the estimated phase aberration $\\{t_{i}\\}$ to the STA RF signals,

$s_{i}(t)=\sum_{j=1}^{128}{P_{ij}{(t-g_{i,j}(x_{f},z_{f})}+t_{i}+t_{j})}.$ (1)

Notice that the phase aberration correction was applied to both $j_{th}$
transmit element and $i_{th}$ receive element. We assumed that each element
has experienced the same phase aberration in both the transmission and
receiving process.

A plot of $s_{i}(t)$ was shown in Fig. 1(a). The focal point was at
$x_{f}=0.96cm$ and $z_{f}=1.63cm$ as signified by the red square in Fig. 1(b).
$s_{i}$ was the signal received at the $i_{th}$ receiving elements after
considering the geometrical delay when a beam focused at $x_{f}$ and $z_{f}$
was transmitted. In this paper depth and time were used interchangeably.
Additionally, we defined the echoes due to the scattering at the $x_{f}$ and
$z_{f}$ as wavefronts. For example, we considered the signals within a window
centered around the depth of 1.63 cm at Fig. 1(a) as wavefronts. When there is
no phase aberration, the echoes from the scattering at $x_{f}$ and $z_{f}$
will make a straight line (see Fig. 1(a) for an example). In reality, there is
always phase aberration, so the wavefront will be distorted like Fig. 7(a).
Our goal was to find the proper delay error that could resolve this
distortion. Here, we introduced a method to estimate the phase aberration
iteratively [13] as follows. We started with an initial guess of phase
aberration $\\{t_{i}^{0}\\}$ equal to zero, then obtain the TBF RF data based
on the current $t_{j}^{k-1}$.

$s_{i}^{k}(t)=\sum_{j=1}^{128}{P_{ij}{(t-g_{i,j}(x_{f},z_{f})}+t_{i}^{k-1}+t_{j}^{k-1})}$
(2)

Here, $s_{i}^{k}(t)$ represented the corrected TBF RF data at the $k_{th}$
iteration. $t_{i}^{k-1}$ was the accumulated phase aberration estimated in the
previous ($k-1$) iterations for the $i_{th}$ receive element. Next, NCC or
f-NCC were applied to $s_{i}^{k}(t)$ to obtain the residual phase aberration
$\\{r_{i}^{k}\\}$ in the $k_{th}$ iteration. See details of finding the
residual phase aberration $\\{r^{k}_{i}\\}$ in the appendix A. Then the
$\\{t_{i}^{k}\\}$ was updated as: $t_{i}^{k}=t_{i}^{k-1}+r_{i}^{k}$ . After
that the updated $\\{t_{i}^{k}\\}$ was applied to Eq. 2 for the next iteration
until $\\{r_{i}^{k}\\}$ converged to zero.

###  II-B Filtered NCC

One of the major shortcomings of NCC techniques is the inadequacy of
conventional NCC method for RF signals with low SNR. As explained in the
introduction section, this requirement can introduce a major problem in STA
data in experimental and in-vivo studies. Therefore it is worthwhile to
investigate a solution.

In addition to the electronic noise, off-axis scattering can interfere with
the wavefront of interest. For example, the oblique wavefront in Fig. 1(a) was
due to the scattering of the off-axis beams (side-lobes or grating lobes) from
the hyper-echoic lesion on the left side of the speckle at $x=0.98$ cm and
$z=1.63cm$ in Fig. 1(b). This off-axis clutter interfered with the estimation
of phase aberration. The inaccurate estimations by NCC illustrated in Fig. 4
are due to this interference. The goal of this study was to develop a phase
aberration estimator in STA imaging in the presence of noise and clutter in RF
data. We investigated solutions that would reduce noise and off-axis signal
while preserving the shape of the wavefront. The aim is to consequently
increase the SNR and Signal-clutter-ratio (SCR) of TBF RF signal, $s_{i}(t)$.
A successful filtering method should not change the temporal wavefront
distortion caused by phase aberration.

In this research, a 2D Fourier Transform (FT) in the aperture and temporal
domain was applied to the 2D TBF RF data, $s_{i}(t)$. This yielded a Fourier
spectrum data (Fig. 2(b)):

$s_{i}(t)\quad\quad-2DFT\rightarrow\quad\quad S_{{k}}(f)$ (3)

In this formulation, $k$ represents the aperture frequency corresponding to
the aperture dimension, and $f$ is the temporal frequency. $S_{k}(f)$ is the
2D Fourier transform of $s_{i}(t)$. It should be noted that after transmit
beamforming, and even in presence of phase aberration, the wavefront should
exhibit a somewhat horizontal pattern and with only minor distortions. In the
spectrum illustrated in Fig. 2(b), horizontal patterns of the aligned
wavefronts correspond to the low frequency components along the aperture
dimension, such as the bright points around the $k=0$ $cm^{-1}$ axis. Off-axis
clutters (oblique lines the in Fig. 1(a)) correspond to higher frequency
components away from $k=0$ $cm^{-1}$ axis. Similarly, noise that lacks the
horizontal pattern (Fig. 2(a)) had higher frequency components across the
Fourier spectrum (Fig. 2(b)). Hence wavefront and noise/clutter can be
differentiated by the differences in their aperture frequencies. By applying a
low-pass filter in the aperture frequency direction $M_{k}$, clutter and noise
will be reduced significantly. In addition, we found that applying a low-pass
filter along the temporal frequency dimension, $N(f)$, could also improve the
quality of TBF signals by reducing the noise in the RF signals. Therefore, we
constructed a two-dimensional oval-shaped filter (shown in Fig. 2(c)) by
multiplying the temporal filter $N(f)$ with the aperture filter $M_{k}$. The
filtered spectrum is:

$S^{\prime}_{{k}}(f)=S_{{k}}(f).M_{k}.N(f).$ (4)

$M_{k}$ and $N(f)$ were constructed from Blackman filters to guarantee a
gradual cutoff in both the aperture and the temporal frequency directions. An
example of the filtered beamformed spectrum is shown in Fig. 2(d). Finally the
filtered TBF signal spectrum $S^{\prime}_{{k}}(f)$ would undergo an inverse
Fourier transform along both the temporal and aperture frequency dimensions.
As depicted in the Fig. 2(e) upon filtering the noise was reduced
significantly in the filtered TBF signals.

To demonstrate the effect of the filter on removing the off-axis clutter
signals, we applied the filter to the TBF signals in Fig. 1(b). The oblique
lines in Fig 1(a) were artifacts generated by the hyper-echoic region that was
in the vicinity of the focal point (Fig. 1(b)). It can be seen in Fig. 3 that
upon filtration the unwanted oblique lines were significantly reduced. It was
shown in Figs. 2 and 3 that the oval filter was capable of removing the off-
axis scattering and noise artifacts in the TBF RF data since they lack a
horizontal pattern. In summary, the oval filters has successfully preserved
the shape of the aberrated wavefront while significantly reducing both the
noise and clutter in RF signals. NCC will be applied to the filtered TBF RF
data that has an improved SNR and SCR (2(e)) to estimate the phase aberration
in the following results section.

##  III Results

###  III-A Simulation Study

Filed II ultrasound simulation [27] software was utilized to generate
synthetic aperture RF signals. A homogeneous medium with a speed of sound of
1540 $m/s$ was simulated. Transmit and receive events were simulated with a
128-element transducer having a pitch of 0.15 $mm$. The simulation and image
reconstruction was performed in Matlab (The MathWorks Inc., Natick, Ma). To
mimic a phase aberration scenario, the phase aberration delay profile (truth
delay) was created by following the recommendation reported by Dahl et al.
[28] by convolving a set of Gaussian random numbers with a Gaussian function.
Phase aberration was applied to both the transmit and receive process of the
STA to simulate the two-way phase aberration in ultrasound propagation in
tissues. In the simulations, this ground truth of the phase aberration profile
was compared with that estimated from the TBF RF signals. To estimate the
aberration delay, NCC algorithm was used, as described by Eq. 6 to 11 in
Appendix A. To generate noisy data, the STA RF signals were degraded by
applying a random noise to individual channels of data to create $RF_{SNR}$ of
-10dB. Noise was filtered to match the bandwidth of the ultrasound signals.

In the simulation study, in the absence of noise, the filtered-NCC method was
evaluated and compared to NCC as illustrated Fig. 4. In the NCC method a
threshold value ($d_{threshold}$) of $\lambda_{0}/4$ was used. This threshold
value is small enough to avoid unwanted large differences between the adjacent
error values. Fig. 4 shows that the results from NCC (red line) deviates
strongly from the ground truth (blue line) for the elements with an index
smaller than 80 (on the left side of the delay profile) while the results of
f-NCC (the yellow line) agrees well with the ground truth. It is noticed that
the off-axis clutter signals (the oblique lines) in Fig. 1(b) intersect with
the wavefront of interest (around the focal depth of 1.63cm) mainly at
elements with index smaller than 80. Therefore, we believe that the clutter
signals caused the deviation of the results (Fig. 4 ) (red line) from the
ground truth (blue line). As illustrated in Fig. 3, the Oval filter in F-NCC
has eliminated the oblique lines. Therefore, f-NCC has estimated the delay
profile accurately.

Later on, the iterative method was applied to the noisy RF signals when focal
point was located at the same depth ($z=1.63cm$) over an speckle area. Fig. 5
shows that NCC estimation of phase aberration delays are different from the
truth as noise was introduced to the RF signal. Ideally, the difference
between the adjacent elements in delay error should be small and the delay
profile should have a smooth pattern. However, in the presence of noise, NCC
fails and can result in a large difference between the adjacent elements in
delay errors. To improve NCC, threshold values of $\lambda_{0}$ and
$\lambda_{0}/4$ were applied to the adjacent element delay errors as explained
in Eq.9 in appendix A. Fig. 5 illustrates that the NCC along with two
different values of threshold could not estimate the error delays accurately.
As shown in Fig. 2, oval filter has removed the superposed noise from the TBF
RF data. Consequently, f-NCC (shown by red line in Fig. 5) were able to
estimate the aberration delay profile accurately. In Fig. 5, a comparison
between truth, NCC with two threshold values and f-NCC delays profiles proves
the success of this filtering method to resolve the limitation of NCC.

The sizes of the oval filter, especially the cutoff frequency in the aperture
direction, are important factors. Various window sizes for the aperture and
temporal filter were tested. Filters with excessively small cutoff frequencies
would over-flatten and alter the overall shape of the wavefront. This is
caused by the removal of the desirable frequencies in the Fourier domains
around $k=0$ $cm^{-1}$ and $f=0$ $MHz$. On the other hand, For large cutoff
windows, filter would allow passage of unwanted frequencies and therefore the
noise and off-axis scattering reduction would not be successful. In this case,
estimated delay profile would be distorted. A decision for the cutoff
frequency was a trade-off between the smoothness of the profile and preserving
the shape of the wavefront. In this study cut-off frequencies of $18cm^{-1}$
and $12MHz$ for oval filters were used.

Due to the periodic nature of Fourier transform, the edge effect was present
in the delay profile estimation, as shown in Fig. 10 in appendix B. To correct
the distortions in the edges of the delay profile an iterative padding
technique was used. See details in appendix B.

Next, Fig. 6 compares the performance of NCC and f-NCC algorithms in TBF data
of RF signals without noise and when $RF_{SNR}=-10dB$. Fig. 7 shows their
corresponding STA images. Fig. 6(a) shows the noiseless, aberrated TBF RF
data. Fig. 6(b) illustrates that NCC has improved the alignment of wavefront
by correcting the phase aberration. Even though in Fig. 4, NCC delay profile
is different from the truth, especially on the left side of the graph, there
is a general agreement in the rest of the graph. Consequently, NCC delay error
compensation has resulted in an improved wavefront in Fig. 6(b). This
observation is compatible with the reconstructed images in Fig. 7(a) and b. In
Fig. 7(a) the smearing of the hyper-echoic lesion, fill-ins into the hypo-
echoic area, and the out-of-focus point reflector at the depth of 2.5cm are
indications of phase aberration. Even though a misalignment can be seen in
Fig. 6(b), NCC has reduced the phase aberration effect significantly for a
noiseless RF in Fig. 7(b). In this image the lesions and point reflectors are
sharper and has better delineation.

Later on, NCC was investigated for a noisy RF ($RF_{SNR}=-10dB$). The TBF data
from the noisy and aberrated data was plotted in Fig. 6(c). The corrected RF
data were presented in Fig. 6(d), showing the poor performance of NCC for
noisy RF data. As illustrated in Fig. 5, NCC is not capable of estimating a
correct phase aberration delay for noisy RF (yellow and purple lines). The
image reconstructed from the noisy and aberrated data in Fig. 6(c) was shown
in Fig. 7(c). There is obvious evidences of phase aberration and noise in this
image. The image reconstructed from the noisy data that had been corrected
with NCC is shown in Fig. 7(d). In this image, the effect of phase aberration
is still present. Even though the focal points are sharper, the hyper-echoic
lesion is smeared and the contrast is not improved.

Fig. 6(e) shows the TBF data after filtration with the designed oval filter.
The noise has been significantly reduced and wavefront is more apparent. The
f-NCC delay profile was estimated from this TBF RF data and as illustrated in
Fig. 5, it is compatible with the truth profile. A comparison between Fig.
6(d) and Fig. 6(f) indicates the advantage of f-NCC versus NCC methods to
flatten the wavefront. After f-NCC delay compensation, the abberated wavefront
shown in Fig. 6(e) was successfully aligned as illustrated in Fig. 6(f). It
should be noted that the corrected TBF RF data plots was achieved after 2
iterations. The corrected image with f-NCC is plotted in Fig. 7(e). Fig. 7
(c), (d) and (e) shows that f-NCC was successful in correcting the phase
aberration effects from the images: 1) image 7(e) has better contrast and
delineates the lesions better. 2) Fig. 7(c) and (d) show blurring in speckles
whereas this issue is significantly improved after f-NCC correction in image
(e). Additionally, the effect of phase aberration correction on the images was
quantitatively compared by measuring the peak signal-to-noise ratio (PSNR) of
the point reflector and the contrast-to-noise-ratio (CNR) values of the hyper-
echoic region. As shown in Table. I values of PSNR and CNR have supported the
above claims. After the f-NCC corrections, CNRhyper and PSNR have been
improved by 24 and 51%, respectively.

###  III-B Experimental Study

####  III-B1 Experiment Setup/Method

The proposed method was investigated on a tissue-mimicking phantom, CIRS
calibrated phantom (Multi-purpose multi-tissue ultrasound phantom Model 040
GSE) with nominal sound-speed of 1540 $\pm$ 20 m/s. A segment with a width of
2.95cm was scanned. To create a noisy RF data, CIRS phantom was topped with a
3.5-cm thick slab of beef, which had a speed of sound of about 1540-1560 m/s
[29]. This calculated speed is similar to the value for the calibrated
phantom. Therefore it was reasonable to assume the phase aberration caused by
the beef is minor. The area of study within the CIRS phantom contained 3 point
reflectors(Fig. 9(a)). Additionally, three hypo-echoic regions, which are
signified with the red arrows, are present within the beef area. Synthetic
aperture raw RF data was acquired with Verosonics Vantage scanner and a P6-4
phased array probe with a central frequency of 3.47 MHz. Sampling frequency
was four times the central frequency. Later on, the RF signals were up-sampled
to 55.5 MHz in software. The probe was fixated by an articulated arm and
placed directly on top of the phantom. The speckle focal point was located in
the CIRS phantom (under the beef), at the imaging depth of 4.1 cm. The ground
truth of the sound speed $c_{t}$ was assumed to be 1540 m/s.

Unlike the simulation study, it is challenging to know the true phase
aberration of an experimental tissue-mimicking phantom. Therefore, in the
experimental study the proposed phase aberration estimation methodology was
studied under a speed error scenario, in which a wrong sound speed of 1640 m/s
was intentionally used to reconstruct the STA images. In the formulation, we
assumed that the phase aberration other than that caused by the speed error
can be neglected. The advantage of a speed error scenario is that the phase
aberration delay profile can be formulated to quantitatively compare NCC and
f-NCC methods. Additionally, in the speed error scenario phase aberration is a
function of depth, as shown in Eq. 5b. This makes an ideal case to investigate
the f-NCC method for STA data. Since NCC can only estimate the relative
shifting between elements, the phase aberration of the $i_{th}$ element was
derived relative to (by subtracting) that of the element right above the focal
point (point C in Fig. 11 of the appendix). The true delay error is formulated
as below:

$\displaystyle\delta t_{i}=\delta t_{c}\cdot(1+cos\theta_{i})$ (5a)
$\displaystyle\delta
t_{c}=l_{i}*(1-cos\theta_{i})*[\frac{1}{c_{r}}-\frac{1}{c_{t}}]$ (5b)

where $i$ represents the receive element index, $\delta t_{i}$ is the phase
aberration error of the $i_{th}$ element relative to that of the element right
above the focal point, $l_{{}_{i}}$ is the travel distance between the focal
point and receive element $i$, and $\theta_{i}$ is the angel between the
travelled path to the receive element $i$ and the norm to the probe, as shown
in Fig. 11 of the appendix. $c_{t}$ represents the true average speed of sound
and $c_{r}$ is the reconstruction sound speed. $\delta t_{c}$ represents the
part of the $i_{th}$ element phase aberration (relative to the element right
above the focal point) caused directly by the sound speed error.

The $cos\theta_{i}$ term in Eq. 5a is an unexpected term. The $cos\theta_{i}$
term is due to the vertical shift of the focal point caused by the speed error
in the image reconstruction. Considering that $\theta_{i}$ usually small, the
$cos\theta_{i}$ term is close to 1. The derivation of the formulation will be
explained in appendix C. The above formula was validated with Field II
simulations and there was a good agreement between the theory and simulations.

Oval filter with the window size of $11cm^{-1}$ in aperture frequency domain
and $17MHz$ in temporal frequency domain was found as the optimum windows in
this experimental study. Ultimately, the experimental RF was corrected over
transmit and receive.

Fig. 8 shows the importance of the filtration to achieve a proper delay
profile estimation. Signals were originated from a speckle area at the depth
of 4.1cm. In this Fig. $c_{t}$ was assumed to be 1540 m/s. The NCC method was
coupled with two values of threshold. The delay profile corresponding to the
threshold value of $\lambda_{0}/4$ exhibits undesirable jumps in the profile.
Nevertheless, it still shows a somewhat similar trend to the true delay errors
in the second half of the delay profile (orange line) indicated by the green
arrow. This suggests the possibility that a smaller threshold value could fix
the problem. However NCC coupled with a smaller threshold ($\lambda_{0}/8$)
shows a substantial difference from the theory. This plots proves the success
of f-NCC estimation when compared to the theory profile (Eq. 5b) with 1540 m/s
as the true speed. The small differences between the theory and f-NCC
estimation might be due to the existence of the spatial speed heterogeneity
within the beef layer. These sound speed variations contribute to the overall
phase aberration delay. As indicated earlier, Eq. 5b only formulates the phase
aberration due to the average speed error and not the spatial speed variation
within the beef. To verify this claim, we collected RF data of the CIRS
phantom in the absence of beef layer and estimated the phase aberration for a
speed error scenario. The f-NCC estimated delay profile agreed well with the
theory.

First, a phase aberration-free image (Fig. 9(a)) was reconstructed with the
correct speed of 1540 m/s. Second, the aberrated image with sound speed error
is illustrated in Fig. 9(b). A comparison between Fig. 9(a) and b shows that
phase aberration degradation is present in the STA image with speed error as:
1) the point reflectors within the phantom are unfocused and 2) the three dark
areas within the beef at around depth of 3cm have fill-ins. The dark areas are
shown by red arrows on Fig. 9(a). Images corrected by f-NCC and NCC With the
focal point at 4.1cm deep are shown in Fig. 9(c) and (d), respectively.

The comparison between these images of Fig. 9 are in agreement with the delay
profiles in Fig. 8. F-NCC was more effective than NCC in improving image
quality because in Fig. 9: 1)image (c) exhibits a more improved spatial
resolution, specially for phantom point reflectors at approximately 4 cm and 5
cm; 2)in image (c), beef structural features at around 3cm (shown by red
arrows in image (a) has better contrast and delineation; 3)the interface
between beef and the CIRS phantom in image (c) is sharper while it has a
broken structure in Fig. 9 (d). 4)speckles in image c has a more uniform and
sharper structure than Fig. 9(d), where speckles have a smaller scale and
broken pattern. This conclusion has been confirmed by image quality metrics.
F-NCC corrections has improved the PSNR by 24% whereas this improvement is 5%
after NCC correction.

It should be pointed out that the phase aberration induced by sound speed
error depends on the position of the focal point or pixel of interest. In the
above correction, we chose the focal point to be a speckle at $x=1.1cm$ and
$z=4.1cm$. Therefore, among the three point reflectors the sharpness of the
nearby point reflector at position (2, 4.1) cm was improved the most.
Alternatively, if we find the phase aberration by focusing on speckles around
other point reflectors, the sharpness of the corresponding reflectors will be
improved the most. In practice, it is more efficient to correct the phase
aberration caused by the speed error by adjusting the average speed in the
image reconstruction, rather than by finding the phase aberration point by
point. Here we use the speed error model to demonstrate the ability of the
proposed method to correct phase aberration pixel by pixel in the image.

##  IV Discussion

In the simulation and experimental study we proved the feasibility of f-NCC to
estimate and correct the phase aberration error. A potential application of
this filtering technique is to de-clutter the TBF RF data in B-mode imaging,
prior to receive beamforming. Consequently, phase aberration in B-mode can be
estimated after TBF RF data is de-cluttered. It should be noted that, f-NCC
can determine phase aberration most accurately around the focal point of the
transmission beam. For speckles that is away from the focal point, the
performance of NCC will be compromised. This is due to the decrease of the
correlation between the TBF RF signals from different elements, which is
caused by the large beam width outside the focal zone [30]. In contrast, STA
imaging will not suffer this limitation because every pixel in STA image is
optimally focused. Nevertheless, the designed filter has the potential to be
used in conjunction with other algorithms that require noise/clutter reduction
as a prior step in B-mode imaging.

Another potential application of Eq. 5b that is under study is to develop an
algorithm to estimate the average speed of a medium. With f-NCC proven to be
an effective estimator of the delay profile $t_{i}$, this equation can be
rearranged to solve for the scalar value $c_{t}$. This problem can easily be
solved and with a low computational load.

Another potential application of the f-NCC algorithm is to develop an
algorithm to map the local speed of a medium. This can be done based on the
fact that, at any speckle point, the phase aberration corresponding to the
array elements is the result of the travel time change that the signal has
experienced along the travelled path. This path is the line connecting the
pixel/speckle and the array elements. The final delay error at every speckle
can be represented as the integration of the sound speed error over the path
[19]. Once we estimate the phase aberration delay of many speckles in an
ultrasound image, we can write down a linear equation set that express the
phase aberration in terms of the local speed error. After that, the standard
methods in solving inverse problems can be applied to obtain the local speed.
We are developing a method along this approach to map the local speed of
tissues with STA.

##  V Conclusion

We demonstrated the efficacy and robustness of an iterative f-NCC method in
STA to estimate the phase aberration at speckle area from single-frame noisy
RF signals in both simulation and experimental study. The 2D filtration of the
TBF data was performed in the aperture and temporal domains to partially
remove off-axis signals and noise. An iterative zero-padding technique was
developed to correct the edge effect caused by filtering. A new formula of the
phase aberration profile was derived and validated when a wrong average sound
speed was used in the image reconstruction. The proposed method can improve
the image quality of STA images with phase aberration. Alternative potential
applications of the method include average speed estimation and local speed
mapping with STA.

##  Appendix A Using Normalized-Cross-correlation to estimate the residual
phase aberration $r_{{i}}$

Normalized-Cross-correlation is a window matching function that quantifies the
relative shift between two similar signals $s_{i}$ & $s_{i+p}$ [8]:

$NCC(\tau)=\frac{\sum_{t=w_{1}}^{w_{2}}s_{{i}}(t)s_{i+p}(t+\tau)}{\sqrt{\sum_{t=w_{1}}^{w_{2}}s_{{i}}^{2}(t)\sum_{t=t_{1}}^{t_{2}}s_{i+p}{{}^{2}(t+\tau)}}}$
(6)

Here $w_{2}$ & $w_{1}$ represent the starting and ending point of the window
in the signals $s_{i}$ & $s_{i+p}$, over which cross correlation is measured.
Additionally, $\tau$ is the time shift between the two signals. Parameter $p$
is the receive element lag between the cross-correlated signals (e.g., $p=1$
indicates the neighboring elements). The shift $\hat{d}_{i,i+p}$, between the
receive elements with lag $p$ is at the peak of the cross-correlation
function:

$\hat{d}_{i,i+p}=\underset{\tau}{\arg\max}(NCC(\tau))$ (7)

In this implementation, signals were up-sampled to 80 MHz and sub-sampling
accuracy in NCC algorithm was achieved with cosine curve fitting. NCC was
performed between a pair of receive elements with lag 1 (adjacent elements) &
lag 2.

In theory, we have

$\hat{d}_{{i+1},{i}}=r_{{i+1}}-r_{{i}},$ (8)

where $r_{{i}}$ and $r_{i+1}$ are the residual delay errors corresponding to
the $i_{th}$ and $(i+1)_{th}$receive elements .

In this paper, we assume that the residual phase aberration profile $r_{{i}}$
is a smooth function of $i$. Therefore, we expect $\hat{d}_{{i}}$ to be small
value compared to the wavelength according to Eq. 8. However, the noise in the
RF signals would cause a drastic increase in some adjacent delay error
calculations and jeopardize an ideally smooth delay profile. A straightforward
solution was to define a threshold for $\hat{d}$:

$\hat{d_{i}}=0\quad if\quad\hat{d_{i}}>d_{threshold}.$ (9)

Eq. 8 can be written into a matrix form as below:

$\quad\quad\quad\quad\quad M\quad\quad\quad\quad\quad\quad\quad\quad
r\quad\quad\quad\quad\quad\quad d$
$\left[\begin{array}[]{cccccc}{1}&{-1}&{0}&{0}&{\cdots}&{0}\\\
{0}&{1}&{-1}&{0}&{\cdots}&{0}\\\ {0}&{0}&{1}&{-1}&{\cdots}&{0}\\\
{\vdots}&{}\hfil&{}\hfil&{}\hfil&{\vdots}\\\
{0}&{0}&{\cdots}&{0}&{1}&{-1}\end{array}\right]\left[\begin{array}[]{c}{r_{1}}\\\
{r_{2}}\\\ {r_{3}}\\\ {\vdots}\\\
{r_{N-1}}\end{array}\right]=\left[\begin{array}[]{c}{d_{1,2}}\\\ {d_{2,3}}\\\
{d_{3,4}}\\\ {\vdots}\\\ {d_{N-1,N}}\end{array}\right]$ (10)

The true arrival time error was estimated by taking the pseudo-inverse of the
above equation:

$\hat{r}=\left(M^{T}M\right)^{-1}M^{T}d$ (11)

The estimation can be improved if M includes both lag 1 and lag 2 shift. That
can be achieved when the corresponding matrices of $M$ and $d$ for both lags
are stacked.

In summary, Eq. 6 was first used to yield the shift between a pair of receive
elements $d_{i},_{i+p}$ with the desired lag $p$. The threshold value was
applied according to equation 9. Then the delay error for each receive element
can be computed by solving Eq. 11. After that, the estimated delay was
compensated to the raw RF data for both the transmit and receive process to
approximately correct the phase aberration. Lastly, the above process was
iterated several times for further improvement. This technique takes
approximately 2-3 iterations to reach convergence and proper correction of the
phase aberration error.

##  Appendix B Padding

In the estimation of the delay profile, a much larger error can be found at
the edge in f-NCC method, as shown by the black arrows at the elements 1-10
and 118-128 in Fig. 10. This is caused by the periodic assumption of signals
when applying Fourier transform. Based on the periodic assumption, signal from
channel 128 will be wrapped round to be adjacent to the signal from channel 1.
Since the signals from these two channels are quite different, there is a
discontinuity in the 2D TBF RF signals. After going through the low-pass
filters, the locations with discontinuity will result in larger distortion
compared to other parts. To accommodate such problems, we used an approach
borrowed from the zero-padding technique in signal processing field [31]. An
iterative padding scheme was proposed. In our implementation of the filtering
process, the 2D TBF RF signals were stored in a matrix $S_{i}(t)$, each column
represented the signal from one channel. In the first iteration, TBF signal
matrix was padded by replicating several times, the signal of the first
channel on the left side of the matrix and the signal of the $128_{th}$
channel on the right side of the matrix. Then this expanded matrix was
filtered as described in the Section II A. After that, f-NCC was applied to
calculate the delay profile. In further iterations, signals were extrapolated
on both sides with delays that have the same slope on the delay profile as
these of the first and last 10 elements. Fig. 10 illustrates that this
iterative padding technique has minimized the edge effect.

##  Appendix C Formulation of Delay Profile under a Speed Error Scenario

To validate the f-NCC methodology, the true phase aberration delay was derived
for a scenario when image was reconstructed with a wrong average speed. The
speed error will result in travel time error in two ways: first, it will
affect the travel time by the wrong speed directly. Second, it will shift the
target vertically in the image, which in turn will cause an additional travel
time error. Initially, for the focal point A (Fig. 11), the travel time
relative to the element right above the focal point was calculated. Later on,
$l_{AC}$ was subtracted so element right above the focal pixel has no delay :

$t_{c}=\frac{l_{AD}-l_{AC}}{c}=\frac{l_{i}}{c}*(1-cos\theta_{i})$ (12)

For the simplicity in this formulation, the path with respect to point $A$ was
shown by $l_{i}$. In the above equation, $\theta_{i}$ is the angel between
travel path with respect to receive element $i$ and norm to the path. The goal
of this calculation was to find the delay error time corresponding to element
$i$, located at point $D$. Under a speed error condition, the delay error is:

$\delta t_{c}=l_{i}*(1-cos\theta_{i})*[\frac{1}{c_{r}}-\frac{1}{c_{t}}]$ (13)

$c_{r}$ represents the reconstruction speed and $c_{t}$ is the true speed.
Additionally, in presence of speed error the focal point A is shifted by
$l_{AB}$ or $\delta l$:

$\delta
l=l_{AB}=(l_{i}-\frac{l_{i}}{c_{t}}.c_{r})*cos\theta_{i}=l_{i}.c_{r}*[\frac{1}{c_{r}}-\frac{1}{c_{t}}]*cos\theta_{i}$
(14)

Time error due to the vertical shift will be:

$\delta t_{d}=[\frac{l_{DB}-l_{BC}}{c_{r}}-\frac{l_{AD}-l_{AC}}{c_{r}}]$ (15)

This equation was rearranged as below :

$\delta t_{d}=[\frac{l_{DB}-l_{AD}}{c_{r}}-\frac{l_{BC}-l_{AC}}{c_{r}}]$ (16)

By substituting for the position shift $\delta l$ this equation was obtained:

$\displaystyle\delta t_{d}=\frac{-\delta lcos(\theta_{i})+\delta
l}{c_{r}}=\frac{\delta l[1-cos(\theta_{i})]}{c_{r}}=$ (17)
$\displaystyle\frac{l_{i}*c_{r}*cos\theta_{i}[1-cos(\theta_{i})]}{c_{r}}*[\frac{1}{c_{r}}-\frac{1}{c_{t}}]=\delta
t_{c}.cos\theta_{i}$

To find the total delay error, speed error $\delta t_{c}$ and position shift
$\delta t_{d}$ delays were added:

$\delta t_{i}=\delta t_{d}+\delta t_{c}=\delta t_{c}.(1+cos\theta_{i})$ (18)

## Acknowledgment

The authors would like to thank the following funding agencies: Natural
Sciences and Engineering Research Council of Canada (NSERC), Canada Foundation
for Innovation (CFI), and Ryerson University. We are grateful to Dr. Michael
Kolios, Dr. Ying Li, Dr. Ping Gong, and Na Zhao for their valuable discussions
and suggestions. We thank Dr. Micheal Kolios for providing access to his
research facility in this study.

## References

  * [1] M. Tabei, T. Mast, and R. Waag, “Simulation of ultrasonic focus aberration and correction through human tissue,” _Journal of the Acoustical Society of America_ , vol. 113, no. 2, pp. 1166–1176, 2003.
  * [2] D. Guenther and W. Walker, “Generalized cystic resolution: A metric for assessing the fundamental limits on beamformer performance,” _IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control_ , vol. 56, no. 1, pp. 77–90, 2009.
  * [3] J. Mann and W. Walker, “A constrained adaptive beamformer for medical ultrasound: Initial results,” vol. 2, 2002, pp. 1807–1810.
  * [4] M. Lediju Bell, J. Dahl, and G. Trahey, “Resolution and brightness characteristics of short-lag spatial coherence (slsc) images,” _IEEE transactions on ultrasonics, ferroelectrics, and frequency control_ , vol. 62, pp. 1265–1276, July 2015.
  * [5] G. Chau, J. Dahl, and R. Lavarello, “Effects of phase aberration and phase aberration correction on the minimum variance beamformer,” _Ultrasonic imaging_ , vol. 40, pp. 16–17, July 2017.
  * [6] M. O’Donnell and S. W. Flax, “Phase-Aberration Correction Using Signals From Point Reflectors and Diffuse Scatterers: Measurements,” _IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control_ , vol. 35, no. 6, pp. 768–774, 1988.
  * [7] S. Krishnan, P.-C. Li, and M. O’Donnell, “Adaptive compensation of phase and magnitude aberrations,” _IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control_ , vol. 43, no. 1, pp. 44–55, 1996.
  * [8] M. O’Donnell and S. W. Flax, “Phase-Aberration Correction Using Signals From Point Reflectors and Diffuse Scatterers: Basic Principales,” _IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control_ , vol. 35, no. 6, pp. 768–774, 1988.
  * [9] D. Liu and R. C. Waag, “Correction of ultrasonic wavefront distortion using backpropagation and a reference waveform method for time-shift compensation.” _The Journal of the Acoustical Society of America_ , vol. 96 2 Pt 1, pp. 649–60, 1994.
  * [10] G. Ng, S. Worrell, P. Freiburger, and G. Trahey, “A comparative evaluation of several algorithms for phase aberration correction,” _IEEE Transactions on Ultrasonics, Ferroelectrics and Frequency Control_ , vol. 41, no. 5, pp. 631–643, 1994. [Online]. Available: http://ieeexplore.ieee.org/lpdocs/epic03/wrapper.htm?arnumber=308498
  * [11] L. Nock, G. E. Trahey, and S. W. Smith, “Phase aberration correction in medical ultrasound using speckle brightness as a quality factor,” _The Journal of the Acoustical Society of America_ , vol. 85, no. 5, pp. 1819–1833, 1989.
  * [12] M. Fink, “Time reversal of ultrasonic fields—part i: Basic principles,” _IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control_ , vol. 39, no. 5, pp. 555–566, 1992.
  * [13] M. Karaman, A. Atalar, H. Koymen, and M. O’Donnell, “A phase aberration correction method for ultrasound imaging.” _IEEE transactions on ultrasonics, ferroelectrics, and frequency control_ , vol. 40, no. 4, pp. 275–282, 1993.
  * [14] B. Byram, K. Dei, J. Tierney, and D. Dumont, “A model and regularization scheme for ultrasonic beamforming clutter reduction,” _IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control_ , vol. 62, no. 11, pp. 1913–1927, 2015.
  * [15] E. Shaswary, J. Tavakkoli, and Y. Xu, “A new algorithm for time-delay estimation in ultrasonic echo signals [Correspondence],” _IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control_ , vol. 62, no. 1, pp. 236–241, 2015.
  * [16] F. Viola and W. F. Walker, “A comparison of the performance of time-delay estimators in medical ultrasound,” _IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control_ , vol. 50, no. 4, pp. 392–401, 2003.
  * [17] B. C. Byram, G. E. Trahey, and J. A. Jensen, “A method for direct localized sound speed estimates using registered virtual detectors,” _Ultrasonic Imaging_ , vol. 34, no. 3, pp. 159–180, 2012.
  * [18] M. E. Anderson and G. E. Trahey, “The direct estimation of sound speed using pulse–echo ultrasound,” _The Journal of the Acoustical Society of America_ , vol. 104, no. 5, pp. 3099–3106, 1998. [Online]. Available: http://asa.scitation.org/doi/10.1121/1.423889
  * [19] M. Jakovljevic, S. Hsieh, R. Ali, G. Chau Loo Kung, D. Hyun, and J. J. Dahl, “Local speed of sound estimation in tissue using pulse-echo ultrasound: Model-based approach,” _The Journal of the Acoustical Society of America_ , vol. 144, no. 1, pp. 254–266, 2018. [Online]. Available: http://asa.scitation.org/doi/10.1121/1.5043402
  * [20] T. Slack, C. Myers, C. Martin, and S. Heymsfield, “The geographic concentration of us adult obesity prevalence and associated social, economic, and environmental factors,” _Obesity (Silver Spring, Md.)_ , vol. 22, March 2014.
  * [21] A.-M. Almeida, H.-P. Cotrim, D.-B. Barbosa _et al._ , “Fatty liver disease in severe obese patients: Diagnostic value of abdominal ultrasound,” _World Journal of Gastroenterology_ , vol. 14, p. 1415, 03 2008.
  * [22] J. Dahl and G. Trahey, “Off-axis scatterer filters for improved aberration measurements,” vol. 2, 11 2003, pp. 343– 347 Vol.1.
  * [23] J. Dahl and T. Feehan, “Direction of arrival filters for improved aberration estimation,” _Ultrasonic imaging_ , vol. 30, pp. 1–20, 02 2008.
  * [24] J. Shin, L. Huang, and J. Yen, “Spatial prediction filtering for medical ultrasound in aberration and random noise,” _IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control_ , vol. PP, pp. 1–1, Oct 2018\.
  * [25] J. Dahl, D. Hyun, Y. Li, M. Jakovljevic, M. Lediju Bell, W. Long, N. Bottenus, V. Kakkad, and G. Trahey, “Coherence beamforming and its applications to the difficult-to-image patient,” Sep 2017, pp. 1–10.
  * [26] P. Gong, M. C. Kolios, and Y. Xu, “Delay-encoded transmission and image reconstruction method in synthetic transmit aperture imaging,” _IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control_ , vol. 62, no. 10, pp. 1745–1756, 2015.
  * [27] J. A. Jensen, S. I. Nikolov, K. L. Gammelmark, and M. H. Pedersen, “Synthetic aperture ultrasound imaging,” _Ultrasonics_ , vol. 44, pp. e5 – e15, 2006\.
  * [28] J. J. Dahl, D. A. Guenther, and G. E. Trahey, “Adaptive Imaging and Spatial Compounding in the Presence of Aberration,” _IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control_ , vol. 52, no. 7, pp. 1131–1144, 2005.
  * [29] S. A. Goss, R. L. Johnston, and F. Dunn, “Comprehensive compilation of empirical ultrasonic properties of mammalian tissues,” _Journal of the Acoustical Society of America_ , vol. 64, no. 2, pp. 423–457, 1978.
  * [30] R. Mallart and M. Fink, “The van Cittert-Zernike theorem in pulse echo measurements,” _Journal of the Acoustical Society of America_ , vol. 90, no. 5, pp. 2718–2727, 1991.
  * [31] A. C. Kak, M. Slaney, and G. Wang, “Principles of computerized tomographic imaging,” _Medical Physics_ , vol. 29, no. 1, pp. 107–107, 2002.

(a)

(b)

Figure 1:  (a) plot of $s_{i}(t)$, transmit beamformed RF signals,
corresponding to the speckle inside the red box in (b). (b) the reconstructed
image in the absence of noise and aberration. Focal depth was at 1.63cm.
a$s_{i}(t)$b$S_{k}(f)$2D
FFTc$M_{k}.N(f)$d$S^{\prime}_{k}(f)$S*Me$s^{\prime}_{i}(t)$iFFT Figure 2:
aperture domain filtering: (a) transmit beamformed signal with a focal point
at speckle point at z=1.63cm. beamforming was performed by all the 128
elements. (b) 2D Fourier transform of the transmit beamformed RF signals. (c)
Low-pass filter along aperture and temporal frequency domains. (d) Filtered RF
in Fourier domain (e) The final result of filtration with the improved and de-
cluttered transmit beamformed data Figure 3:  Transmit beamformed plots after
filtration. The off-axis scattering is significalty reduced Figure 4: A
comparison between the truth, NCC and f-NCC delay estimation in absence of
noise. In the NCC method threshold was equal to $\lambda_{0}/4$. Focal point
was at speckle area at 1.63cm depth. Focal point was in vicinity of a hyper-
echoic region.  Figure 5: A comparison between the truth, NCC, and f-NCC delay
estimation. Independent noise was added to the ultrasound signal to achieve a
${SNR}$ of -10 dB. NCC method was coupled with thresholds values of
$\lambda_{0}$ and $\lambda_{0}/4$ for adjacent element values. Focal point was
on a speckle region at the depth of 1.63cm  No NoiseAberratedaCorrectedb-10db
NoisecdFilteredef Figure 6:  Plots of Transmit beamformed aberrated RF signals
(a) with no noise, (c) with added noise of -10dB SNR and (e) after filtration
of the noisy signals in (c). (b, d) and (f) show their corresponding plots
after phase aberration correction with NCC. Focal point was at a speckle point
at z=1.63 cm.  No NoiseAberratedaCorrectedbSNR=-10dB cdNCCef-NCC Figure 7:
STA images reconstructed from : (a) the noiseless aberrated data, (b) the NCC-
corrected noiseless data, (c) noisy RF signals with aberration, (d) the NCC-
corrected noisy data, and (e) the f-NCC-corrected noisy data. Focal point was
located at speckle point at z=1.63 cm Figure 8:  A comparison between the
theory and the estimated phase aberration. NCC coupled with thresholds of
$\lambda_{0}/4$ and $\lambda_{0}/8$ and f-NCC were included. Theory delay
profile was calculated using Eq. 5b with respect to focal point at 4.1 cm.

(a)

(b)

(c)

(d)

Figure 9: Experimental STA images reconstructed from RF signals with (a) a correct speed of 1540 m/s and (b) a wrong speed of 1650 m/s. STA images reconstructed from the RF signals corrected with (c) f-NCC and (d) NCC methods. A beef top layer was added to create a noisy RF data. The focal point to estimate the phase aberration (signified by the red square) was at 4.1 cm. Figure 10: A comparison between the truth and estimated delay profiles to show the effectiveness of padding to correct the edge effect. In this graph blue line represents the Truth delay profile. Red and yellow lines compare the corrected profiles with oval filter with and without padding. Figure 11: The illustration of truth travel delay calculation. Points D shows the transducer element and point A is the focal point position. C is the transducer element right about the focal point A  TABLE I: Quantitative comparison of aberrated and NCC and f-NCC corrected images of simulated phantom | PSNR [dB] | CNR hyper [dB]
---|---|---
Aberrated | 5.4 | 1.54
NCC | 5.8 | 1.2
f-NCC | 8.2 | 1.92


# Higher-Order Topological Insulator on a Martini Lattice and Its Square Root
Descendant

Daiki Matsumoto1 Tomonari Mizoguchi2 and Yasuhiro Hatsugai2 1Graduate School
of Pure and Applied Sciences1Graduate School of Pure and Applied Sciences
University of Tsukuba University of Tsukuba Tsukuba Tsukuba Ibaraki 305-8571
Ibaraki 305-8571 Japan 2Department of Physics Japan 2Department of Physics
University of Tsukuba University of Tsukuba Tsukuba Tsukuba Ibaraki 305-8571
Ibaraki 305-8571 Japan Japan

###### Abstract

The notion of square-root topological insulators has been recently generalized
to higher-order topological insulators. In two-dimensional square-root higher-
order topological insulators, the emergence of in-gap corner states is
inherited from the squared Hamiltonian, which hosts higher-order topology. In
this paper, we first propose that the martini lattice model serves as a
concrete example of a higher-order topological insulator. We further propose a
square-root higher-order topological insulator based on the martini lattice
model. Specifically, we propose that the honeycomb model with two-site
decoration, whose squared Hamiltonian consists of two martini lattice models,
realizes square-root higher-order topological insulators. We show, for both
the martini lattice and the decorated honeycomb model, that in-gap corner
states appear at finite energies and that non-trivial bulk $\mathbb{Z}_{3}$
topological invariant protects them.

##  1 Introduction

The topological phase of matter is one of the hot topics in condensed matter
physics. Topological insulators(TIs) are the typical example of the
topological phase. [1, 2, 3, 4, 5, 6]. They have robust gapless states at
$(d-1)$-dimensional boundaries of samples due to non-trivial topological
indices($d$ is the dimension of bulk). This relation between bulk topological
indices and boundary states is known as bulk-boundary correspondence [7, 8],
which is one of the most characteristic nature of the topological phase.

Nowadays, various kinds of TIs have been actively explored. Higher-order
topological insulators (HOTIs) are one of the examples of such novel
topological phases [9, 10, 11, 12, 13, 14]. HOTIs in $d$ dimensions have
$(d-n)$-dimensional topologically-protected boundary states (with $n\geq 2$)
rather than $(d-1)$-dimensional topological boundary states as conventional
TIs. For instance, a second-order TI in two dimensions has $0-$dimensional
boundary states, i.e, the corner states[15, 16, 17, 18].

Figure 1: (Color online)Schematic figures of (a) the martini lattice model,
and (b) the honeycomb lattice model with two-site decorations. Here
$\bm{a}_{1}$ and $\bm{a}_{2}$ represent the lattice vectors:
$\bm{a}_{1}=\left(\frac{1}{2},\frac{\sqrt{3}}{2}\right)$ and
$\bm{a}_{2}=\left(-\frac{1}{2},\frac{\sqrt{3}}{2}\right)$. (c) An illustration
of the relation between the original Hamiltonian and its square. The squared
Hamiltonian of the decorated honeycomb lattice model is the direct sum of two
martini lattice models with sublattice-dependent on-site potentials.

Square-root TIs are another example of a novel topological phase [19].
Historically, the concept of square-root plays an essential role in various
aspects of physics [20, 21, 22, 23, 24, 25]. In the context of the topological
phase, square-root TIs stand for the TIs that inherit their non-trivial
topological nature from the square of the Hamiltonian, which is referred to as
a parent Hamiltonian. Consequently, a Hamiltonian for a square-root TI is
referred to as a child Hamiltonian. Since its discovery, experimental
realizations of square-root TIs have actively been pursued, mainly in
artificial topological materials [26, 27, 28, 29, 30, 31, 32] such as photonic
crystals, phononic crystals, and electric circuits. On the theoretical side,
the generation of square-root descendant is applied to various topological
phases [33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46]. One such
example is a generalization to HOTIs, which is dubbed as square-root HOTIs
[33]. Since the square-root HOTIs are a novel phase and only a few examples
are known so far, it is desirable to search for additional concrete examples
to verify their ubiquity.

In this paper, we first propose a tight-binding model which hosts a HOTI, that
is, the martini lattice model [47, 48, 49] shown in Fig. 1(a). Then, based on
this, we also propose a realization of square-root HOTI in the honeycomb
lattice model with two-site decoration [50, 51], which we simply call the
decorated honeycomb lattice model in what follows, shown in Fig. 1(b). The
squared Hamiltonian of the decorated honeycomb model corresponds to the direct
sum of two martini lattices as schematically shown in Fig. 1(c). In other
words, the martini lattice model is the parent of the decorated honeycomb
lattice model. Thus, the square-root HOTI is realized in the decorated
honeycomb lattice model. Remarkably, we find that changing the corner
termination changes the nature of the topological corner modes, which
originates from the fact that the corner termination specifies the
corresponding parent Hamiltonian from which the corner modes are inherited.

The rest of this paper is organized as follows. In Sect. 2, we introduce the
martini lattice model and the decorated honeycomb lattice model and clarify
the key symmetries. Then, we point out their square-root relation with a
special focus on the inheritance of the energy dispersion and Bloch wave
functions. We also show the dispersion relations obtained by numerical
calculation for each model to seek topological phase transition points. In
Sect. 3, we elucidate the higher-order topology of the martini and the
decorated honeycomb models, by demonstrating the existence of in-gap corner
states in the finite system under the open boundary conditions. We also
numerically calculate $\mathbb{Z}_{3}$ topological indices related to
polarization for each model to confirm the bulk-corner correspondence. It is
revealed that the decorated honeycomb model realizes the square-root HOTI
inherited from the martini model. In Sect. 4, we present a summary of this
paper.

##  2 Model

###  2.1 Martini model

As a parent Hamiltonian of the following discussions, we consider a tight-
binding Hamiltonian on a martini lattice. This lattice structure has been
considered in several previous works [52, 48, 49]. It has four sites per unit
cell, as shown in Fig. 1(a). The tight-binding Hamiltonian reads

$H^{(\bm{M})}=\sum_{\bm{k}}\bm{C}_{\bm{k}}^{\dagger}H^{(\bm{M})}_{\bm{k}}\bm{C}_{\bm{k}},$
(1)

where
$\bm{C}_{\bm{k}}=(C_{\bm{k},(\rm{i})},C_{\bm{k},(\rm{ii})},C_{\bm{k},(\rm{iii})},C_{\bm{k},(\rm{iv})})^{T}$
are the annihilation operators, and

$H^{(\bm{M})}_{{\bm{k}}}=\left(\begin{array}[]{cccc}V&t_{2}&t_{2}e^{i\bm{k}\cdot\bm{a}_{1}}&t_{2}e^{i\bm{k}\cdot\bm{a}_{2}}\\\
t_{2}&0&t_{1}&t_{1}\\\ t_{2}e^{-i\bm{k}\cdot\bm{a}_{1}}&t_{1}&0&t_{1}\\\
t_{2}e^{-i\bm{k}\cdot\bm{a}_{2}}&t_{1}&t_{1}&0\end{array}\right),$ (2)

is the Hamiltonian matrix. For the definition of $\bm{a}_{1}$ and
$\bm{a}_{2}$, see Fig. 1(a). This model includes three parameters: nearest-
neighbor hoppings with two different parameters, namely, $t_{1}$ for intra-
triangles and $t_{2}$ for inter-triangles, and the difference of on-site
potentials between the black and white sites, $V$. We note, for later
discussions, that the Hamiltonian has $C_{3}$ symmetry centered at the center
of the red triangle consisting of sublattice (i), (ii), and (iii) in a unit
cell: $H^{(\bm{M})}_{\bm{k}}$ satisfies

$H^{(\bm{M})}_{C_{3}\bm{k}}=U^{(\bm{M})}_{\bm{k}}H^{(\bm{M})}_{\bm{k}}(U^{(\bm{M})}_{\bm{k}})^{\dagger},$
(3)

where
$C_{3}\bm{k}=(-\frac{1}{2}k_{x}+\frac{\sqrt{3}}{2}k_{y},-\frac{\sqrt{3}}{2}k_{x}-\frac{1}{2}k_{y})$
and

$U^{(\bm{M})}_{\bm{k}}=\left(\begin{array}[]{cccc}e^{-i\bm{k}\cdot\bm{a}_{2}}&0&0&0\\\
0&0&0&1\\\ 0&1&0&0\\\ 0&0&1&0\end{array}\right).$ (4)

The Hamiltonian has also $C_{3}$ symmetry centered at sublattice(i):

$\displaystyle H^{(\bm{M})}_{C_{3}\bm{k}}$
$\displaystyle=U^{(\bm{dM})}_{\bm{k}}H^{(\bm{M})}_{\bm{k}}(U^{(\bm{dM})}_{\bm{k}})^{\dagger},$
(5)

where

$\displaystyle U^{(\bm{dM})}_{\bm{k}}$
$\displaystyle=e^{-i\bm{k}\cdot\bm{a}_{2}}(U^{(\bm{M})}_{\bm{k}})^{*}.$ (6)

We depict the band structure for several sets of hopping parameters with
fixing $V=0$ in Fig. 2. We see that a flat band appears at $E=-t_{1}$
regardless of parameters, which originates from the fact that the martini
lattice belongs to a class of lattices called partial line graphs[47]. When
$t_{1}/t_{2}<0$, the third band is a flat band [Fig. 2 (a)-(d)], whereas when
$t_{1}/t_{2}>0$, the second band is a flat band [Fig. 2 (e)-(h)]. In general,
we see two gaps in the band structure and we focus on these gaps when
discussing the topological phases. The band touching between the first and the
second bands occurs for $t_{1}/t_{2}=-1/\sqrt{2}$, and that between the third
and the fourth bands occurs for $t_{1}/t_{2}=1/\sqrt{2}$; they occur at K
point. Additionally, the triple band touching among the second, the third, and
the fourth bands occur for $t_{1}/t_{2}=-1$, and that among the first, the
second, and the third bands occur for $t_{1}/t_{2}=1$; they occur at $\Gamma$
point.

Figure 2: (Color online)Band structure of the martini lattice model with
$t_{2}=1$,$V=0$ and (a) $t_{1}/t_{2}=-1.5$, (b) $t_{1}/t_{2}=-1$, (c)
$t_{1}/t_{2}=-1/\sqrt{2}$, (d) $t_{1}/t_{2}=-0.3$, (e) $t_{1}/t_{2}=0.3$, (f)
$t_{1}/t_{2}=1/\sqrt{2}$, (g) $t_{1}/t_{2}=1$, (h) $t_{1}/t_{2}=1.5$.
$\Gamma=(0,0),{{\rm{K}}}=\left(\frac{4\pi}{3},0\right)$ and
${\rm{M}}=\left(\pi,\frac{\pi}{\sqrt{3}}\right)$ are high-symmetry points in
the Brillouin zone.

###  2.2 Decorated honeycomb model

Next, we introduce a tight-binding model on a decorated honeycomb lattice.
This lattice structure has been considered theoretically [50, 51], and is also
relevant to some solid-state materials [53, 54, 55, 56]. The Hamiltonian is

$H^{(\bm{DH})}=\sum_{\bm{k}}\bm{C}_{\bm{k}}^{\dagger}H^{(\bm{DH})}_{\bm{k}}\bm{C}_{\bm{k}},$
(7)

where
$\bm{C}_{\bm{k}}=(C_{\bm{k},\bullet(\rm{I})},C_{\bm{k},\bullet(\rm{II})},C_{\bm{k},\bullet(\rm{III})},C_{\bm{k},\bullet(\rm{IV})},\\\
C_{\bm{k},\circ(\rm{i})},C_{\bm{k},\circ(\rm{ii})},C_{\bm{k},\circ(\rm{iii})},C_{\bm{k},\circ(\rm{iv})})^{T}$
are the annihilation operators, and

$H^{(\bm{DH})}_{\bm{k}}=\left(\begin{array}[]{cc}\mathcal{O}_{4,4}&\Phi_{\bm{k}}^{\dagger}\\\
\Phi_{\bm{k}}&\mathcal{O}_{4,4}\end{array}\right),$ (8)

is the Hamiltonian matrix. In Eq. (8), we have introduced

$\Phi_{\bm{k}}=\left(\begin{array}[]{cccc}0&t_{2}&t_{2}&t_{2}\\\
t_{1}&t_{3}&0&0\\\ t_{1}&0&t_{3}e^{i\bm{k}\cdot\bm{a}_{1}}&0\\\
t_{1}&0&0&t_{3}e^{i\bm{k}\cdot\bm{a}_{2}}\end{array}\right),$ (9)

and $n\times m$ zero matrix, $\mathcal{O}_{n,m}$. Note that this model
includes only nearest-neighbor hoppings with three different parameters
$t_{1},t_{2}$ and $t_{3}$. The vectors $\bm{a}_{1}$ and $\bm{a}_{2}$ are
defined in the same way as for a martini lattice model.

As for the symmetries of this Hamiltonian, we note that the Hamiltonian has
chiral symmetry:

$\gamma H^{(\bm{DH})}_{\bm{k}}\gamma=-H^{(\bm{DH})}_{\bm{k}},$ (10)

where

$\gamma=\left(\begin{array}[]{cc}I_{4}&\mathcal{O}_{4,4}\\\
\mathcal{O}_{4,4}&-I_{4}\end{array}\right).$ (11)

Here $I_{n}$ is a $n$ dimensional identity matrix. The Hamiltonian also has
$C_{3}$ symmetry centered at sublattice $\circ$(i):

$H^{(\bm{DH})}_{C_{3}\bm{k}}=U^{(\bm{DH})}_{\bm{k}}H^{(\bm{DH})}_{\bm{k}}(U^{(\bm{DH})}_{\bm{k}})^{\dagger},$
(12)

where

$U^{(\bm{DH})}_{\bm{k}}=\left(\begin{array}[]{cc}U^{(\bm{M})}_{\bm{k}}&\mathcal{O}_{4,4}\\\
\mathcal{O}_{4,4}&U^{(\bm{dM})}_{\bm{k}}\end{array}\right).$ (13)

Figure 3: (Color online)Band structures of decorated honeycomb lattice with
$t_{1}=1,t_{3}=\sqrt{\frac{3}{5}}$ and (a) $t_{2}=0.4$, (b) $t_{2}=0.5$, (c)
$t_{2}=0.8$, (d) $t_{2}=1$, (e) $t_{2}=1.2$.
$\Gamma=(0,0),{{\rm{K}}}=\left(\frac{4\pi}{3},0\right)$ and
${\rm{M}}=\left(\pi,\frac{\pi}{\sqrt{3}}\right)$ are high-symmetry points in
the Brillouin zone.

In Fig. 3, we show the band structures of the decorated honeycomb model. Here,
and also in Sect. 3, we fix the parameters $t_{1}$ and $t_{3}$ as
$\left(t_{1},t_{3}\right)=\left(1,\sqrt{3/5}\right)$, and retain $t_{2}$ as
the only variable parameter, for the purpose of concreteness. We see that two
flat bands appear for all the panels of Fig. 3. In fact, the flat bands appear
at $E=\pm t_{3}$ regardless of the choice of parameters[51]. In Fig. 3(b), we
see that triple band touching among the second, the third, and the fourth
bands and that among the fifth, the sixth, and the seventh bands occur at
$\Gamma$ point with $t_{2}=1/2$. In fact, it is analytically shown that the
triple band touching at $\Gamma$ point occurs for
$(1/t_{1}^{2}+1/t_{2}^{2})t_{3}^{2}=3$ [51]. Besides $\Gamma$ point, the band
touchings between the first and the second bands as well as the seventh and
the eighth bands occur at K point for $|t_{1}|=|t_{2}|$ [Figs. 3(d)]. For
$|t_{1}|\neq|t_{2}|$ they are gapped, and we focus on these gaps when
discussing the square-root HOTIs, in the decorated honeycomb model.

###  2.3 Square-root relation

We now show the square-root relation between the decorated honeycomb lattice
model and the martini lattice model. A simple understanding of this relation
can be obtained from the real-space viewpoint. Since the Hamiltonian includes
only NN hoppings which connect only white sites with black sites, a particle
on a black (white) site can moves to the nearest black (white) sites or comes
back to the original site by operating the Hamiltonian twice. The former paths
network black (white) sites that form a martini (downward martini) lattice,
while the latter corresponds to the on-site potentials. Turning to the square-
root relation in the reciprocal space, i.e., that for the Bloch Hamiltonian,
we see that the square of the decorated honeycomb lattice model becomes

$[H^{(\bm{DH})}_{\bm{k}}]^{2}=\left(\begin{array}[]{cc}h_{\bm{k}}^{(\bm{M})}&\mathcal{O}_{4,4}\\\
\mathcal{O}_{4,4}&(h_{\bm{k}}^{(\bm{M})})^{\dagger}\end{array}\right),$ (14)

where

$\displaystyle h_{\bm{k}}^{(\bm{M})}$ $\displaystyle=$
$\displaystyle\Phi_{\bm{k}}^{\dagger}\Phi_{\bm{k}}$ (15) $\displaystyle=$
$\displaystyle\left(\begin{array}[]{cccc}3t_{1}^{2}&t_{1}t_{3}&t_{1}t_{3}e^{i\bm{k}\cdot\bm{a}_{1}}&t_{1}t_{3}e^{i\bm{k}\cdot\bm{a}_{2}}\\\
t_{1}t_{3}&t_{2}^{2}+t_{3}^{2}&t_{2}^{2}&t_{2}^{2}\\\
t_{1}t_{3}e^{-i\bm{k}\cdot\bm{a}_{1}}&t_{2}^{2}&t_{2}^{2}+t_{3}^{2}&t_{2}^{2}\\\
t_{1}t_{3}e^{-i\bm{k}\cdot\bm{a}_{2}}&t_{2}^{2}&t_{2}^{2}&t_{2}^{2}+t_{3}^{2}\end{array}\right),$
(20)

is equivalent to the Hamiltonian of a martini lattice of Eq. (2), with the
nearest-neighbor hoppings being $t_{1}t_{3}$ and $t_{2}^{2}$, and the on-site
potential being $3t_{1}^{2}-\left(t_{2}^{2}+t_{3}^{2}\right)$, up to the
constant energy shift $t_{2}^{2}+t_{3}^{2}$. Further, for the lower-right part
of Eq. (14), we have

$\displaystyle(h_{\bm{k}}^{(\bm{M})})^{\dagger}$ $\displaystyle=$
$\displaystyle\Phi_{\bm{k}}\Phi_{\bm{k}}^{\dagger}$ (21) $\displaystyle=$
$\displaystyle\left(\begin{array}[]{cccc}3t_{2}^{2}&t_{2}t_{3}&t_{2}t_{3}e^{-i\bm{k}\cdot\bm{a}_{1}}&t_{2}t_{3}e^{-i\bm{k}\cdot\bm{a}_{2}}\\\
t_{2}t_{3}&t_{1}^{2}+t_{3}^{2}&t_{1}^{2}&t_{1}^{2}\\\
t_{2}t_{3}e^{i\bm{k}\cdot\bm{a}_{1}}&t_{1}^{2}&t_{1}^{2}+t_{3}^{2}&t_{1}^{2}\\\
t_{2}t_{3}e^{i\bm{k}\cdot\bm{a}_{2}}&t_{1}^{2}&t_{1}^{2}&t_{1}^{2}+t_{3}^{2}\end{array}\right).$
(26)

Importantly, $(h_{\bm{k}}^{(\bm{M})})^{\dagger}$ is also the Hamiltonian of
the martini lattice with the upside-down orientation to that of
$h_{\bm{k}}^{(\bm{M})}$ [see Fig. 1(c)]. The corresponding tight-binding
parameters are given as follows: the nearest-neighbor hoppings are
$t_{2}t_{3}$ and $t_{1}^{2}$, and the on-site potential is
$3t_{2}^{2}-\left(t_{1}^{2}+t_{3}^{2}\right)$, and the constant energy shift
is $t_{1}^{2}+t_{3}^{2}$. We call the lattice for
$(h_{\bm{k}}^{(\bm{M})})^{\dagger}$ as a downward martini lattice and write
$(h^{(\bm{M})}_{\bm{k}})^{\dagger}$ as $h^{(\bm{dM})}_{\bm{k}}$ henceforth.
Note that the rotation operator $U^{(\bm{dM})}_{\bm{k}}$ in Eq. (13) acts on
the subspace of the downward martini lattice. Based on this relation, we
elucidate that the energy dispersion and the eigenstates of the decorated
honeycomb lattice can be expressed by using those of squared Hamiltonian [33].
Let $\ket{u_{\bm{k},i}^{(\bm{M})}}$ be the normalized eigenstate of
$h_{\bm{k}}^{(\bm{M})}$ with the eigenenergy
$E^{(\bm{M})}_{{\bm{k}},i}$($i=1,2,3,4$). Note that we label the eigenstate
such that $E^{({\bm{M}})}_{{\bm{k}},1}\leq\cdots\leq
E^{({\bm{M}})}_{{\bm{k}},4}$ holds. Namely, $\ket{u_{\bm{k},i}^{(\bm{M})}}$
satisfies

$h_{\bm{k}}^{(\bm{M})}\ket{u_{\bm{k},i}^{(\bm{M})}}=\Phi_{\bm{k}}^{\dagger}\Phi_{\bm{k}}\ket{u_{\bm{k},i}^{(\bm{M})}}=E^{(\bm{M})}_{{\bm{k}},i}\ket{u_{\bm{k},i}^{(\bm{M})}}.$
(27)

Then, multiplying $\Phi_{\bm{k}}$ on Eq. (27) from left, we get

$\displaystyle\Phi_{\bm{k}}\Phi_{\bm{k}}^{\dagger}(\Phi_{\bm{k}}\ket{u_{\bm{k},i}^{(\bm{M})}})$
$\displaystyle=$ $\displaystyle
h^{(\bm{dM})}_{\bm{k}}\Phi_{\bm{k}}\ket{u_{\bm{k},i}^{(\bm{M})}}$ (28)
$\displaystyle=$ $\displaystyle
E^{(\bm{M})}_{{\bm{k}},i}\Phi_{\bm{k}}\ket{u_{\bm{k},i}^{(\bm{M})}}.$

This relation means that $h_{\bm{k}}^{(\bm{dM})}$ has the same eigenenergies
as $h_{\bm{k}}^{(\bm{M})}$, and that the corresponding eigenstate is

$\ket{u_{\bm{k},i}^{(\bm{dM})}}=\Phi_{\bm{k}}\ket{u_{\bm{k},i}^{(\bm{M})}}.$
(29)

Using these relations, we are ready to construct the eigenstates of
$H^{({\bm{DH}})}_{{\bm{k}}}$. To be concrete, for the non-zero energy state,
the positive energy eigenstate of $H^{({\bm{DH}})}_{{\bm{k}}}$ is given as

$\ket{u_{\bm{k},i+4}}=\frac{1}{N_{\bm{k},i}}\left(\begin{array}[]{c}\sqrt{E_{\bm{k},i}^{(\bm{M})}}\ket{u_{\bm{k},i}^{(\bm{M})}}\\\
\ket{u_{\bm{k},i}^{(\bm{dM})}}\end{array}\right),$ (30)

with the eigenenergy

$\varepsilon_{\bm{k},i+4}=\sqrt{E^{(\bm{M})}_{\bm{k},i}}.$ (31)

Note that $N_{\bm{k},i}=\sqrt{2E_{\bm{k},i}^{(\bm{M})}}$ is the normalization
constant. One can explicitly check this as

$\displaystyle H^{(\bm{DH})}_{\bm{k}}\ket{u_{\bm{k},i+4}}$ $\displaystyle=$
$\displaystyle\frac{1}{N_{\bm{k},i}}\left(\begin{array}[]{c}\Phi^{\dagger}_{\bm{k}}\ket{u_{\bm{k},i}^{(\bm{dM})}}\\\
\sqrt{E_{\bm{k},i}^{(\bm{M})}}\Phi_{\bm{k}}\ket{u_{\bm{k},i}^{(\bm{M})}}\end{array}\right)$
(34) $\displaystyle=$
$\displaystyle\frac{1}{N_{\bm{k},i}}\left(\begin{array}[]{c}\Phi^{\dagger}_{\bm{k}}\Phi_{\bm{k}}\ket{u_{\bm{k},i}^{(\bm{M})}}\\\
\sqrt{E_{\bm{k},i}^{(\bm{M})}}\Phi_{\bm{k}}\ket{u_{\bm{k},i}^{(\bm{M})}}\end{array}\right)$
(37) $\displaystyle=$
$\displaystyle\frac{1}{N_{\bm{k},i}}\left(\begin{array}[]{c}h^{(\bm{M})}_{\bm{k}}\ket{u_{\bm{k},i}^{(\bm{M})}}\\\
\sqrt{E_{\bm{k},i}^{(\bm{M})}}\ket{u_{\bm{k},i}^{(\bm{dM})}}\end{array}\right)$
(40) $\displaystyle=$
$\displaystyle\frac{1}{N_{\bm{k},i}}\left(\begin{array}[]{c}E_{\bm{k},i}^{(\bm{M})}\ket{u_{\bm{k},i}^{(\bm{M})}}\\\
\sqrt{E_{\bm{k},i}^{(\bm{M})}}\ket{u_{\bm{k},i}^{(\bm{dM})}}\end{array}\right)$
(43) $\displaystyle=$
$\displaystyle\sqrt{E_{\bm{k},i}^{(\bm{M})}}\ket{u_{\bm{k},i+4}},$ (44)

where we use Eq. (29) on the second and third lines. We note that, due to the
chiral symmetry of Eq. (10), the negative-energy eigenstate is given as
$\gamma\ket{u_{\bm{k},i+4}}$, and its eigenenergy is
$-\sqrt{E_{\bm{k},i}^{(\bm{M})}}$. We also note that for the zero-energy state
which occurs at K and $\rm{K}^{\prime}$ points with $i=1$, the above
construction of eigenstate can not be applied because
$E^{({\bm{M}})}_{{\rm{K}},1}=0$. The zero-energy modes are doubly degenerate
and they can be chosen as (for K point)

$\displaystyle\ket{u_{{\rm{K}},4}}=\left(\begin{array}[]{c}\ket{u_{{\rm{K}},1}^{(\bm{M})}}\\\
\bm{0}\end{array}\right),\ket{u_{{\rm{K}},5}}=\left(\begin{array}[]{c}\bm{0}\\\
\ket{u_{{\rm{K}},1}^{(\bm{dM})}}\end{array}\right).$ (49)

Here $\bm{0}$ is a four-component zero vector.

##  3 Topological Corner State

In this section, we analyze the higher-order topology of the martini and the
decorated honeycomb lattice model, putting emphasis on the corner-shape
dependence as well as the consequences of the square-root relation.

###  3.1 Martini lattice

We first argue the higher-order topological phase of the martini model. We
first examine the existence of the corner states, paying particular attention
to the corner-shape dependence of the boundary modes. To this end, we study
two finite samples formed into triangle shape with different corner
termination under the open boundary conditions (OBC), shown in Fig. 4(a1) and
(a2). We call the finite sample shown in Fig. 4(a1) $S_{M}$ and that in 4(a2)
$S_{dM}$. The energy spectrums of $S_{M}$ and $S_{dM}$ as a function of
$t_{1}/t_{2}$ are shown in Fig. 4(b1) and (b2). In Fig. 4(b1), we see that
there appear two in-gap states for $-1<t_{1}/t_{2}<1/2$ (with positive energy)
and for $-1/2<t_{1}/t_{2}<1$ (with negative energy). On the other hand, in
Fig. 4(b2), we see that there appear two in-gap states for
$t_{1}/t_{2}<-1,t_{1}/t_{2}>2$ (with positive energy) and for
$t_{1}/t_{2}<-2,t_{1}/t_{2}>1$ (with negative energy). As we have seen in Fig.
2, $t_{1}/t_{2}=-1$ and $t_{1}/t_{2}=1/\sqrt{2}$ ($t_{1}/t_{2}=-1/\sqrt{2}$
and $t_{1}/t_{2}=1$) are the gap-closing points for the lower (upper) gap,
indicating the topological phase transition at these points. In Fig. 4(c1) and
(c2), we plot the energy spectrums of $S_{M}$ and $S_{dM}$ at
$(t_{1},t_{2})=(0.2,1)$ and $(t_{1},t_{2})=(3,1)$. Green and red ellipses are
the in-gap states with negative and positive energies, respectively. We see
that each in-gap state has three-fold degeneracy, which is because the sample
hosts three corners.

To further reveal the real-space properties of these in-gap states, we plot
the real-space probability density distribution of these in-gap states in Fig.
4(d1) and (d2) for negative energy modes and Fig. 4(e1) and (e2) for positive
energy modes. Note that we take the average over the three degenerate states.
We see that the in-gap states have large amplitude at the corner and hence are
corner states.

Figure 4: (Color online)(a1-2) Finite systems under OBC, consisting of 84
black sites and 36 white sites(a1) and 84 black sites and 28 white sites(a2).
(b1-2) The energy spectrums. The horizontal axis is $t_{1}/t_{2}$.The dashed
line and the solid line correspond to the spectrum of periodic boundary
conditions (PBC) and OBC respectively. The green and red lines correspond to
the in-gap states with negative and positive energies, respectively. (c1) The
energy spectrum of $S_{M}$ for $t_{1}=0.2,t_{2}=1$. (c2) The energy spectrum
of $S_{dM}$ for $t_{1}=3,t_{2}=1$. The in-gap corner states are encircled by
green (lower) and red (upper) ellipses. The probability density distribution
of (d1-2) the lower in-gap states and (e1-2) the upper in-gap states on
$S_{M}$ and $S_{dM}$. The radii of gray circles represent the probability
density. (f) The schematic picture of adiabatic limit of $t_{1}=0$. A
tetramer, a trimer, and a dimer are encircled. (g) The energy levels of the
tetramer, the trimer, and the dimer. (h) The schematic picture of a corner of
$S_{dM}$ when $t_{1}\gg t_{2}$.

It is worth noting that the corner states appear at finite energy. This is in
contrast to the breathing kagome model where the corner state has zero
energy[13]. This can be understood by the following discussion of adiabatic
limit[13]. Namely, we focus on the case of $t_{1}=0$ which can be connected to
HOTI phase without closing the band gap. When $t_{1}=0$, $S_{M}$ consists of
the tetramers in the bulk, the trimers at the edges, and the three dimers at
the corners as shown in Fig. 4(f), whereas $S_{dM}$ consists of the tetramers.
The tetramers have the energy $\pm\sqrt{3}t_{2}$ and two-fold zero energy, the
trimers have the energy $\pm\sqrt{2}t_{2}$ and zero energy, and the dimers
have the energy $\pm t_{2}$ [Fig. 4(g)]. As the energy levels of corner dimers
are separated from those of tetramers and trimers, these dimers serve as in-
gap corner states. This holds even for finite $t_{1}$ as far as the band gap
remains open, which accounts for the emergence of finite-energy corner states
in $S_{M}$. The appearance of the corner modes of $S_{dM}$ is understood as
follows. When $t_{1}\gg t_{2}$, at a corner of the sample, two sites are not
connected to any of the red ($t_{1}$) bonds. These two sites again form a
dimer with the hopping $t_{2}$, resulting in the corner states with $E=\pm
t_{2}$. See Fig. 4(h) for the schematic figure of the above argument.

Next, we discuss the topological origin of the in-gap corner states. Here, we
use the $\mathbb{Z}_{3}$ topological invariant $p_{n}$, which is related to
polarization[57, 13, 33]. For $C_{3}$-symmetric system, it takes the form:

$2\pi p_{n}=\arg{\xi_{n}({\rm{K}})}\ (\rm{mod}\ 2\pi),$ (50)

where

$\xi_{n}(\bm{k})=\bra{u_{\bm{k},n}}U_{\bm{k}}\ket{u_{\bm{k},n}}.$ (51)

Here $\ket{u_{{\bm{k}},n}}$ is an eigenstate of Hamiltonian $H_{\bm{k}}$ with
band index $n$ and $U_{\bm{k}}$ is a $C_{3}$ operator that satisfies
$H_{C_{3}{\bm{k}}}=U_{{\bm{k}}}H_{{\bm{k}}}U^{\dagger}_{{\bm{k}}}$. At K
point, $H_{\bm{k}}$ and $U_{\bm{k}}$ are commutable since K point is a
$C_{3}$-invariant momentum. This results in a $\mathbb{Z}_{3}$ quantization of
$p_{n}$, $p_{n}=l/3$ with $l=0,1,2$. More generally, when the target bands
have degeneracy points and the isolated band is not well-defined on such
points, we can define the $p$ for a set of $N$ bands
$\mathcal{N}=\\{N_{1},N_{2},\cdots,N_{N}\\}$ as

$2\pi p_{\mathcal{N}}=\arg{[\det{\Xi({\rm{K}})}]}\ (\rm{mod}\ 2\pi),$ (52)

where $\Xi_{i,j}(\bm{k})$ is the $N\times N$ matrix whose $i,j$ component is
given as

$\Xi_{i,j}(\bm{k})=\bra{u_{\bm{k},i}}U_{\bm{k}}\ket{u_{\bm{k},j}}\
(i,j\in\mathcal{N}).$ (53)

To characterize the lower in-gap corner states, we numerically calculate
$p_{1}$ for $-2<t_{1}/t_{2}<1$ and $p_{\\{1,2\\}}$ for $1<t_{1}/t_{2}<2$.
Similarly, to characterize the higher in-gap corner states, we numerically
calculate $p_{\\{3,4\\}}$ for $-2<t_{1}/t_{2}<1$ and $p_{4}$ for
$1<t_{1}/t_{2}<2$. We plot the numerical result in Fig. 5. In Fig. 5 (a1) and
(b1), we use $U^{(\bm{M})}_{\bm{k}}$ as the $C_{3}$ operator, on the other
hand, in Fig. 5(a2) and (b2), we use $U^{(\bm{dM})}_{\bm{k}}$ as the $C_{3}$
operator.

Figure 5: (Color online)Numerical results of the polarization p of Eq. (50)
and (52) as a function of $t_{1}/t_{2}$. (a1-2) Blue circles are for the first
band. Orange triangles are for the first and second bands. (b1-2) Blue circles
are for the third fourth band. Orange triangles are for the fourth band.

Figure 6: (Color online)(a1-2) The energy spectrums of $S_{M}$ and $S_{dM}$
with $V=0.5$. The dashed line and the solid line correspond to the spectrum of
PBC and OBC respectively. The green and red lines correspond to the in-gap
states with negative and positive energies, respectively. The probability
density distribution of (b1-2) lower in-gap states and (c1-2) upper in-gap
states in (a1-2).

We see that $p$ is equal to $1/3$ for the parameters where the in-gap states
appear, namely, $-1/2<t_{1}/t_{2}<1$ ($t_{1}/t_{2}<-2,t_{1}/t_{2}>1$) for
lower in-gap states and $-1<t_{1}/t_{2}<1/2$ ($t_{1}/t_{2}<-1,t_{1}/t_{2}>2$)
for higher in-gap states of $S_{M}$ ($S_{dM}$). The value of $p$ jumps at
$t_{1}/t_{2}=-1/\sqrt{2}$ and $t_{1}/t_{2}=1$ for $p_{1}$ and $p_{\\{1,2\\}}$
and at $t_{1}/t_{2}=-1$ and $t_{1}/t_{2}=1/\sqrt{2}$ for $p_{4}$ and
$p_{\\{3,4\\}}$, where the band-gap closes and topological transition occurs.
These results establish the bulk-corner correspondence of the HOTI in the
martini model.

So far, we deal with the case of $V=0$, but we point out that the in-gap
corner states appear even when $V\neq 0$. Our numerical results regarding the
corner states are shown in Fig. 6. This result is essential for considering
the square-root relation to the decorated honeycomb model, as we have shown in
Eqs. (15) and (21), as well as we discuss in the next subsection.

###  3.2 Decorated honeycomb lattice

Now we proceed to study the topological nature of the decorated honeycomb
model. To this end, we study finite samples. We study two finite samples with
different corner terminations[Fig. 7(a) and 8(a)], to investigate how the
corner states are inherited from the parents. In the following, we call the
finite sample of Fig. 7(a) [Fig. 8(a)] $S_{\textit{1}}$ [$S_{\textit{2}}$].

Before going to the results, we remark on the reason why we adopt these two
corner terminations. To explain this, it is important to point out the parent-
child relation under the OBC. As is the case of the PBC, the squared
Hamiltonian for the decorated honeycomb model is expressed as a direct sum of
the upward martini (composed of black sites) and the downward martini
(composed of white sites) models [Fig. 1(c)]. However, unlike the PBC case,
the sites on the edges and corners acquire additional on-site potentials after
the squared operation, because those sites have fewer coordination numbers in
the child lattice. Then, the corner termination for $S_{\textit{1}}$
($S_{\textit{2}}$) is chosen such that the no black (white) site on the edges
and corners has a coordination-number imbalance while letting such imbalance
for the white (black) sites. In other words, the corner termination is chosen
such that one of the two parent Hamiltonians becomes free from the un-
uniformity at the boundaries. Consequently, the “natural” parent for a given
corner termination is the other block that does not have un-uniformity. To be
concrete, for $S_{\textit{1}}$ ($S_{\textit{2}}$), the upward martini model of
$S_{M}$ (the upward martini model of $S_{dM}$) is the natural parent.

Based on the above insight, we argue the corner modes and the topological
characterization. We first argue the results in Fig. 7, i.e., those for
$S_{1}$. We see that there appear two types of in-gap states represented by
the green line and red line for $|t_{2}|<1/2$ and $|t_{2}|<1$. As shown in
Fig. 3, $|t_{2}|=1$ and $|t_{2}|=1/2$ are the gap-closing points at $\Gamma$
point and K point, indicating the topological phase transitions occur at these
points. In Fig. 7 (c), we show the energy spectrum for
$(t_{1},t_{2},t_{3})=(1,0.2,\sqrt{3/5})$. Green and red ellipses indicate the
lower and upper in-gap states. We see that each in-gap state has three-fold
degeneracy, which is again because the sample hosts three corners. We plot the
real-space probability density distribution of the in-gap states encircled by
a green ellipse in Fig. 7 (d). We also plot that of the in-gap states
encircled by a red ellipse in Fig. 7 (e).

Figure 7: (Color online)(a) Finite system under OBC, consisting of 120 black
sites and 136 white sites of the decorated honeycomb model. (b) The energy
spectrum with $t_{1}=1,t_{3}=\sqrt{3/5}$. The horizontal axis is $t_{2}$. The
dashed line and the solid line correspond to the spectrum of PBC and OBC
respectively. The green and red lines correspond to the in-gap states with
negative and positive energies, respectively. (c) The energy spectrum of lower
half energies for $t_{1}=1,t_{2}=0.2,t_{3}=\sqrt{3/5}$. The in-gap corner
states are encircled by red and green ellipses. The probability density
distribution of (d) the in-gap states encircled by the red ellipse ($S_{1}-1$)
and (e) the in-gap states encircled by the green ellipse ($S_{1}-2$) in (c).
The radii of gray circles represent the probability density.

We again take the average over the three degenerate states. We see that, in
both Figs. 7(d) and 7(e), the in-gap states have large amplitude at the corner
and hence are corner states. Meanwhile, there is a difference between Figs.
7(d) and 7(e) with respect to sublattice dependence. Namely, the corner state
of the bottom left corner of Fig. 7(d) has a large amplitude at sublattice
(iii), on the other hand, that of the bottom left corner of Fig. 7(e) has a
large amplitude at sublattice (iii) and sublattice (I). We call the in-gap
states shown in Fig. 7(d) [Fig. 7(e)] $S_{\textit{1}}$-1 [$S_{\textit{1}}$-2].

Figure 8: (Color online)(a) Finite system under OBC, consisting of 120 black
sites and 112 white sites of the decorated honeycomb model. (b) The energy
spectrum with $t_{1}=1,t_{3}=\sqrt{3/5}$. The horizontal axis is $t_{2}$. The
dashed line and the solid line correspond to the spectrum of PBC and OBC
respectively. The green and red lines correspond to the in-gap states with
negative and positive energies, respectively. (c1) The energy spectrum for
$t_{1}=1,t_{2}=0.3,t_{3}=\sqrt{3/5}$. The in-gap corner states are encircled
by green ellipses.(c2) The energy spectrum for
$t_{1}=1,t_{2}=1.0,t_{3}=\sqrt{3/5}$. The in-gap corner states are encircled
by red ellipses. The probability density distribution of (d1) the lower in-gap
states in (c1) and (d2) the lower in-gap states in (c2). The radii of gray
circles represent the probability density.

Figure 9: (Color online)Numerical results of the polarization $p$ of Eq. (50)
and (52) as a function of $t_{2}$. The blue circles in these pictures are for
the corresponding parent Hamiltonian. (a) The $p$ characterize
$S_{\textit{1}}$-1. Orange squares are for the first band of the decorated
honeycomb model. (b) The $p$ characterize $S_{\textit{1}}$-2. Orange squares
are for the fourth and fifth bands of the decorated honeycomb model. Green
rhombuses are for the third, fourth, fifth, and sixth bands of the decorated
honeycomb model. (c) The $p$ characterize $S_{\textit{2}}$-1. Orange squares
are for the first band of the decorated honeycomb model. (d) The $p$
characterize $S_{\textit{2}}$-2. Orange squares are for the first, second, and
third bands of the decorated honeycomb model. Green rhombuses are for the
first and second bands of the decorated honeycomb model.

Similarly, we also see the corner states for the sample $S_{\textit{2}}$ , as
shown in Fig. 8. Interestingly, in contrast to $S_{\textit{1}}$, we see that
the corner states, the red corner states ($S_{\textit{2}}$-1) and the green
corner states ($S_{\textit{2}}$-2) do not coexist in the same parameter
region. Rather, $S_{\textit{2}}$-1 appears for $|t_{2}|\lesssim 1$ while
$S_{\textit{2}}$-2 appears for $|t_{2}|\gtrsim 1$ [Fig. 8(b)], which can be
accounted for by the corresponding result for the parent Hamiltonian.

We discuss the topological origin of the in-gap corner states. Here, we use
$p$ of Eqs. (50) and (52) as a topological invariant. We characterize the in-
gap states $S_{\textit{1}}$-1 and $S_{\textit{1}}$-2 by the polarization of
the decorated honeycomb model and that of the upward martini model. To be
concrete, to characterize the in-gap state $S_{\textit{1}}$-1, we numerically
calculate $p_{1}$ and corresponding $p$ of upward martini model
$p^{(\bm{M})}$. To characterize the in-gap state $S_{\textit{1}}$-2, we also
numerically calculated $p_{\\{4,5\\}}$ for $|t_{2}|<1/2,p_{\\{3,4,5,6\\}}$ for
$|t_{2}|>1/2$ and corresponding $p$ of the upward martini model
$p^{(\bm{M})}$. Similarly, we characterize the in-gap states
$S_{\textit{2}}$-1 and $S_{\textit{2}}$-2 by the polarization of the decorated
honeycomb model and that of the downward martini model. To be concrete, to
characterize the in-gap state $S_{\textit{2}}$-1, we numerically calculate
$p_{1}$ and corresponding $p$ of the downward martini model $p^{(\bm{dM})}$.
To characterize the in-gap state $S_{\textit{2}}$-2, we also numerically
calculated $p_{\\{1,2\\}}$ for $|t_{2}|<1/2,p_{\\{1,2,3\\}}$ for $|t_{2}|>1/2$
and corresponding $p$ of of downward martini model $p^{(\bm{dM})}$. We plot
the numerical result in Fig. 9. We see that the $p$ is equal to $1/3$ for the
parameters where the in-gap states appear, namely, $|t_{2}|<1$ for
$S_{\textit{1}}$-1 and $S_{\textit{2}}$-1, and $|t_{2}|>1/2$ for
$S_{\textit{1}}$-2 and $S_{\textit{2}}$-2. The value of $p$ jumps at
$|t_{2}|=1$ for $p_{1}$ and $|t_{2}|=1/2$ for $p_{\\{1,2\\}}$,
$p_{\\{1,2,3\\}}$, $p_{\\{4,5\\}}$, $p_{\\{3,4,5,6\\}}$, where the band-gap
closes and topological transition occurs. the $p$s for the parents and the
children are identical, except for the region $|t_{2}|>1/2$ in Fig. 9(b). In
Appendix, we explain the origin of this behavior. We emphasize that this does
not mean the breakdown of the inheritance of topology from the parent to the
child.

The results in Fig. 9 establish the bulk-corner correspondence of the HOTI in
the decorated honeycomb model. The results also indicate that the emergence of
the corner states on the decorated honeycomb model is inherited from the
martini model, which means the realization of square-root HOTI on the
decorated honeycomb model.

##  4 Summary

In this paper,we have proposed a concrete example of the HOTI and its square-
root. Namely, we have shown that the conventional HOTI is realized in the
martini lattice model, and the square-root HOTI is realized in the decorated
honeycomb lattice with two sites on each edge of the hexagon. The emergence of
the corner states and their topological protection by the bulk
$\mathbb{Z}_{3}$ invariant have been confirmed.

Remarkably, we find that the corner states in the decorated honeycomb model
are sensitive to the corner termination. From the viewpoint of the square-root
topology, this can be accounted for by the fact that the un-uniformity around
the boundary appears in one of the blocks of the squared Hamiltonian:
Consequently, the “natural” parent for a given corner termination is the other
block which do not have un-uniformity.

On the experimental realization, we believe that both the martini lattice and
the decorated honeycomb lattice are feasible in artificial topological
systems, such as photonic crystals, phononic crystals, and electric circuits.
In particular, for the decorated honeycomb model, the energy of corner state
can be switched by varying $t_{2}$ [Fig. 8(b)], which might serve as
intriguing property for a corner-mode-based-engineering.

## Acknowledgement

This work is supported by JSPS KAKENHI, Grands No. JP17H06138, No. JP20K14371
(T.M.) and JST-CREST JPMJCR19T1(Y.H.).

## Appendix: Relation of the polarization between parents and child

We point out the generic relation of the polarization between the parents and
the child. Note that the case of $E^{(\bm{M})}_{{\bm{k}},i}=0$ is an
exception, which occurs at K point with $i=1$.

We first argue the non-zero energy modes, in particular, for the positive-
energy modes of the decorated honeycomb lattice model. A similar argument
holds for the negative energy sector. The positive-energy eigenvector of Eq.
(30) and the block diagonal form of the $C_{3}$ operator of Eq. (13) lead to

$\displaystyle U^{(\bm{DH})}_{\rm{K}}\ket{u_{{\rm{K}},i+4}}$
$\displaystyle=\frac{1}{\sqrt{2E_{{\rm{K}},i}^{(\bm{M})}}}\left(\begin{array}[]{c}\sqrt{E_{{\rm{K}},i}^{(\bm{M})}}U^{(\bm{M})}_{\rm{K}}\ket{u_{{\rm{K}},i}^{(\bm{M})}}\\\
U^{(\bm{dM})}_{\rm{K}}\ket{u_{{\rm{K}},i}^{(\bm{dM})}}\end{array}\right)$ (56)
$\displaystyle=\frac{1}{\sqrt{2E_{{\rm{K}},i}^{(\bm{M})}}}\left(\begin{array}[]{c}\xi^{(\bm{M})}_{i}({\rm{K}})\sqrt{E_{{\rm{K}},i}^{(\bm{M})}}\ket{u_{{\rm{K}},i}^{(\bm{M})}}\\\
\xi^{(\bm{dM})}_{i}({\rm{K}})\ket{u_{{\rm{K}},i}^{(\bm{dM})}}\end{array}\right),$
(A.1)

where $\xi_{i}^{(\bm{M})}({\rm{K}})$ and $\xi_{i}^{({\bm{dM}})}({\rm{K}})$ are
the eigenvalues of $U^{(\bm{M})}_{{\rm{K}}}$ and $U^{{\bm{dM}}}_{{\rm{K}}}$
for $\ket{u^{({\bm{M}})}_{{\rm{K}},i+4}}$ and
$\ket{u^{({\bm{dM})}}_{{\rm{K}},i+4}}$, respectively. Since
$\ket{u_{{\rm{K}},i+4}}$ is the eigenstate of $U_{{\rm{K}}}^{(\bm{DH})}$, the
eigenvalue $\xi_{i+4}({\rm{K}})$ has to be identical to
$\xi_{i}^{({\bm{M}})}({\rm{K}})$ and $\xi_{i}^{({\bm{dM}})}({\rm{K}})$.
Consequently, for $p$, the relation

$\displaystyle p_{n}=p^{(\bm{M})}=p^{(\bm{dM})},$ (A.2)

holds. Indeed, we see that this relation holds in Fig. 9, except for Fig.
9(b).

In fact, in Fig. 9(b), we consider the polarization including the case of
$E^{({\bm{M}})}_{{\rm{K}},i}=0$, corresponding to the fourth and the fifth
bands for $H_{{\rm{K}}}^{({\bm{DH}})}$. As we have mentioned in Sect. 2.3, the
zero-energy eigenstate of $H^{({\bm{DH}})}_{{\rm{K}}}$ can be chosen as Eq.
(49). Both of these two states are the eigenstates of
$U_{{\rm{K}}}^{({\bm{DH}})}$, but their eigenvalues are, respectively,
$\xi_{1}^{({\bm{M}})}({\rm{K}})$, and $\xi_{1}^{({\bm{dM}})}({\rm{K}})$, which
are not necessarily identical in this case. Then, unlike the case of finite-
energy states, we can argue the parent-child correspondence is only for a set
of these two bands. Specifically, $p_{\\{4,5\\}}$ satisfies

$\displaystyle p_{\\{4,5\\}}=p^{({\bm{M}})}_{1}+p^{({\bm{dM}})}_{1}.$ (A.3)

In Fig. A.1, we plot the $t_{2}$ dependence of
$p^{({\bm{dM}})},p^{({\bm{M}})}$, and $p_{\\{4,5\\}}$ /
$p_{\\{3,4,5,6\\}}$corresponding to the choice of bands of Fig. 9(b). We see
that the relation of Eq. (A.3) indeed holds.

Figure A.1: (Color online)Numerical results of the polarization $p$ of Eq.
(50) and (52) as a function of $t_{2}$. Blue circles are for the corresponding
parent Hamiltonian. Orange squares are for the fourth and fifth bands of the
decorated honeycomb lattice model. Green rhombuses are for the third, fourth,
fifth, and sixth bands of the decorated honeycomb lattice model.

## References

  * [1] F. D. M. Haldane: Phys. Rev. Lett. 61 (1988) 2015.
  * [2] C. L. Kane and E. J. Mele: Phys. Rev. Lett. 95 (2005) 146802.
  * [3] C. L. Kane and E. J. Mele: Phys. Rev. Lett. 95 (2005) 226801.
  * [4] B. A. Bernevig, T. L. Hughes, and S.-C. Zhang: Science 314 (2006) 1757\.
  * [5] M. Z. Hasan and C. L. Kane: Rev. Mod. Phys. 82 (2010) 3045.
  * [6] X.-L. Qi and S.-C. Zhang: Rev. Mod. Phys. 83 (2011) 1057.
  * [7] Y. Hatsugai: Phys. Rev. Lett. 71 (1993) 3697.
  * [8] Y. Hatsugai: Phys. Rev. B 48 (1993) 11851.
  * [9] W. A. Benalcazar, B. A. Bernevig, and T. L. Hughes: Science 357 (2017) 61.
  * [10] W. A. Benalcazar, B. A. Bernevig, and T. L. Hughes: Phys. Rev. B 96 (2017) 245115.
  * [11] K. Hashimoto, X. Wu, and T. Kimura: Phys. Rev. B 95 (2017) 165443.
  * [12] F. Schindler, A. M. Cook, M. G. Vergniory, Z. Wang, S. S. P. Parkin, B. A. Bernevig, and T. Neupert: Science Advances 4 (2018) eaat0346.
  * [13] M. Ezawa: Phys. Rev. Lett. 120 (2018) 026801.
  * [14] S. Hayashi: Communications in Mathematical Physics 364 (2018) 343.
  * [15] H. Araki, T. Mizoguchi, and Y. Hatsugai: Phys. Rev. B 99 (2019) 085406\.
  * [16] Y. Takane: Journal of the Physical Society of Japan 88 (2019) 094712\.
  * [17] H. Watanabe and S. Ono: Phys. Rev. B 102 (2020) 165120.
  * [18] R. Takahashi, T. Zhang, and S. Murakami: Phys. Rev. B 103 (2021) 205123\.
  * [19] J. Arkinstall, M. H. Teimourpour, L. Feng, R. El-Ganainy, and H. Schomerus: Phys. Rev. B 95 (2017) 165109.
  * [20] P. A. M. Dirac and R. H. Fowler: Proceedings of the Royal Society of London. Series A, Containing Papers of a Mathematical and Physical Character 117 (1928) 610.
  * [21] C. L. Kane and T. C. Lubensky: Nature Physics 10 (2014) 39.
  * [22] J. Attig and S. Trebst: Phys. Rev. B 96 (2017) 085145.
  * [23] J. Attig, K. Roychowdhury, M. J. Lawler, and S. Trebst: Phys. Rev. Research 1 (2019) 032047.
  * [24] G. G. Naumis, L. A. Navarro-Labastida, E. Aguilar-Méndez, and A. Espinosa-Champo: Phys. Rev. B 103 (2021) 245418.
  * [25] L. A. Navarro-Labastida, A. Espinosa-Champo, E. Aguilar-Mendez, and G. G. Naumis: Phys. Rev. B 105 (2022) 115434.
  * [26] M. Kremer, I. Petrides, E. Meyer, M. Heinrich, O. Zilberberg, and A. Szameit: Nature Communications 11 (2020) 907.
  * [27] L. Song, H. Yang, Y. Cao, and P. Yan: Nano Letters 20 (2020) 7566.
  * [28] M. Yan, X. Huang, L. Luo, J. Lu, W. Deng, and Z. Liu: Phys. Rev. B 102 (2020) 180102.
  * [29] H. Wu, G. Wei, Z. Liu, and J.-J. Xiao: Opt. Lett. 46 (2021) 4256.
  * [30] W. Yan, D. Song, S. Xia, J. Xie, L. Tang, J. Xu, and Z. Chen: ACS Photonics 8 (2021) 3308.
  * [31] J. Kang, T. Liu, M. Yan, D. Yang, X. Huang, R. Wei, J. Qui, G. Dong, Z. Yang, and F. Nori: arXiv:2109.00879 (2021).
  * [32] Z.-G. Geng, Y.-G. Peng, H. Lv, Z. Xiong, Z. Chen, and X.-F. Zhu: Journal of Physics: Condensed Matter 34 (2021) 104001.
  * [33] T. Mizoguchi, Y. Kuno, and Y. Hatsugai: Phys. Rev. A 102 (2020) 033527 [Erratum 104 (2021) 029906].
  * [34] M. Ezawa: Phys. Rev. Research 2 (2020) 033397.
  * [35] T. Mizoguchi, T. Yoshida, and Y. Hatsugai: Phys. Rev. B 103 (2021) 045136\.
  * [36] A. M. Marques, L. Madail, and R. G. Dias: Phys. Rev. B 103 (2021) 235425\.
  * [37] R. G. Dias and A. M. Marques: Phys. Rev. B 103 (2021) 245112.
  * [38] T. Yoshida, T. Mizoguchi, Y. Kuno, and Y. Hatsugai: Phys. Rev. B 103 (2021) 235130.
  * [39] A. M. Marques and R. G. Dias: Phys. Rev. B 104 (2021) 165410.
  * [40] L. Song, H. Yang, Y. Cao, and P. Yan: Nature Communications 13 (2022) 5601.
  * [41] W. Cheng, X. Zhang, M.-H. Lu, and Y.-F. Chen: Phys. Rev. B 105 (2022) 094103.
  * [42] R.-L. Zhang, Q.-P. Wu, M.-R. Liu, X.-B. Xiao, and Z.-F. Liu: Annalen der Physik n/a 2100497\.
  * [43] K. Roychowdhury, J. Attig, S. Trebst, and M. J. Lawler: arXiv:2207.09475 (2022).
  * [44] R. W. Bomantara: Phys. Rev. B 106 (2022) L060305.
  * [45] W. Deng, T. Chen, and X. Zhang: Phys. Rev. Research 4 (2022) 033109\.
  * [46] L. Zhou, R. W. Bomantara, and S. Wu: SciPost Phys. 13 (2022) 015.
  * [47] S. Miyahara, K. Kubo, H. Ono, Y. Shimomura, and N. Furukawa: Journal of the Physical Society of Japan 74 (2005) 1918.
  * [48] C. R. Scullard: Phys. Rev. E 73 (2006) 016107.
  * [49] P. A. McClarty, M. Haque, A. Sen, and J. Richter: Phys. Rev. B 102 (2020) 224303.
  * [50] C. Barreteau, F. Ducastelle, and T. Mallah: Journal of Physics: Condensed Matter 29 (2017) 465302.
  * [51] T. Mizoguchi, H. Katsura, I. Maruyama, and Y. Hatsugai: Phys. Rev. B 104 (2021) 035155.
  * [52] K. Kubo, C. Hotta, S. Miyahara, and N. Furukawa: Physica B: Condensed Matter 378-380 (2006) 273.  Proceedings of the International Conference on Strongly Correlated Electron Systems.
  * [53] R. H. Baughman, H. Eckhardt, and M. Kertesz: The Journal of Chemical Physics 87 (1987) 6687.
  * [54] R. Longuinhos, E. A. Moujaes, S. S. Alexandre, and R. W. Nunes: Chemistry of Materials 26 (2014) 3701.
  * [55] Z. Li, M. Smeu, A. Rives, V. Maraval, R. Chauvin, M. A. Ratner, and E. Borguet: Nature Communications 6 (2015) 6321.
  * [56] J. M. Lee, C. Geng, J. W. Park, M. Oshikawa, S.-S. Lee, H. W. Yeom, and G. Y. Cho: Phys. Rev. Lett. 124 (2020) 137002.
  * [57] C. Fang, M. J. Gilbert, and B. A. Bernevig: Phys. Rev. B 86 (2012) 115112\.