[type=editor, auid=000,bioid=1, orcid=0009-0000-1982-7110] [*]

[*]

1]organization=UPES, addressline=Bidholi road, city=Dehradun, postcode=248001,
country=India, ,*Equal Contribution

[style=chinese] [*]

[]

2]organization=Ulster, addressline=Northland Road, city=Londonderry,
postcode=BT1 6DN., country=United Kingdom

3]organization=IIIT-Delhi, addressline=Okhla Road, postcode=695571,
state=Delhi, country=India

#  Modality-Order Matters! A Novel Hierarchical Feature Fusion Method for
CoSAm: A Code-Switched Autism Corpus

Mohd Mujtaba Akhtar<EMAIL_ADDRESS>[  Girish  Muskaan Singh  [
Orchid Chetia Phukan  [

###### Abstract

Autism Spectrum Disorder (ASD) is a complex neuro-developmental challenge,
presenting a spectrum of difficulties in social interaction, communication,
and the expression of repetitive behaviors in different situations. This
increasing prevalence underscores the importance of ASD as a major public
health concern and the need for comprehensive research initiatives to advance
our understanding of the disorder and its early detection methods. This study
introduces a novel hierarchical feature fusion method aimed at enhancing the
early detection of ASD in children through the analysis of code-switched
speech (English and Hindi). Employing advanced audio processing techniques,
the research integrates acoustic, paralinguistic, and linguistic information
using Transformer Encoders. This innovative fusion strategy is designed to
improve classification robustness and accuracy, crucial for early and precise
ASD identification. The methodology involves collecting a code-switched speech
corpus, CoSAm, from children diagnosed with ASD and a matched control group.
The dataset comprises 61 voice recordings from 30 children diagnosed with ASD
and 31 from neurotypical children, aged between 3 and 13 years, resulting in a
total of 159.75 minutes of voice recordings. The feature analysis focuses on
MFCCs and extensive statistical attributes to capture speech pattern
variability and complexity. The study employs a range of classical machine
learning models, including SVM, RF, KNN, and so on and deep learning
techniques such as CNN, RNN, and Transformer. Significantly, the research
highlights the importance of modality order in feature fusion, demonstrating
enhanced diagnostic capabilities when specific sequences are used. The best
model performance is achieved using a hierarchical fusion technique with an
accuracy of 98.75% using a combination of acoustic and linguistic features
first, followed by paralinguistic features in a hierarchical manner. These
results underscore the potential of sophisticated signal processing techniques
in the early detection of neurological conditions and contribute valuable
insights into speech processing in bilingual or multilingual individuals with
ASD.``

###### keywords:

Autism Spectrum Disorder (ASD) Hierarchical Feature Fusion Transformer
Encoders Paralinguistic Features Speech Pattern Biomarkers Advanced Signal
Processing

The proposed method is a unique hierarchical fusion strategy that utilizes
Transformer Encoders to successfully integrate acoustic, paralinguistic, and
linguistic information. This methodology aims to improve the identification of
Autism Spectrum Disorder (ASD) in children.

We utilize a novel curated code-switched (English and Hindi) speech corpus and
the proposed feature fusion model achieved remarkable results in detecting ASD
among children in the given corpus exhibiting high classification accuracy and
robustness.

Emphasized the method’s capacity to facilitate earlier and more precise
identification of ASD, a critical factor in developing timely and efficacious
intervention approaches.

Insights into the speech patterns of bilingual or multilingual individuals
with ASD, a relatively unexplored area in speech processing are provided that
will be of great value to future research.

##  1 Introduction

Autism Spectrum Disorder (ASD) stands as a complex neurodevelopmental
challenge, presenting a spectrum of difficulties in social interaction,
communication, and the expression of repetitive behaviors (American
Psychiatric Association [APA], 2013). The term "spectrum" underscores the
varied nature of symptoms experienced by individuals with ASD, reflecting
unique strengths and challenges for each affected individual. ASD McCann and
Peppé (2003); Diehl et al. (2009); Kjelgaard and Tager-Flusberg (2001)
typically reveals itself in early childhood, often surfacing before the age of
three. Early signs may include a lack of responsiveness to social cues,
limited eye contact, delayed speech development, and the manifestation of
repetitive movements or behaviors. The global prevalence of ASD is on the
rise, with the Centers for Disease Control and Prevention (CDC) reporting that
approximately 1 out of 36 children in the United States is identified with
ASD. This escalating prevalence underscores the significance of ASD as a major
public health concern and the need for comprehensive research initiatives to
advance our understanding of the disorder and its early detection methods. In
this context, our research embarks on a detailed exploration of speech
features to contribute to the growing body of knowledge to enhance the timely
and accurate identification of ASD in children. This paper navigates through
dataset collection, preprocessing, and exploratory analysis, setting the stage
for the subsequent application of machine learning models. The synthesized
insights aim to contribute to the growing body of knowledge in the domain of
ASD detection, providing a foundation for more accurate and timely
interventions for affected individuals. In the pursuit of a comprehensive
understanding of ASD detection in children, this research employs advanced
speech processing techniques, focusing on the extraction and analysis of
MFCCs, prosodic features, and statistical attributes. Comprising a total of 39
features, MFCC features are computed using the OpenSMILE toolkit. This
structured storage facilitates subsequent data analysis.

Furthermore, statistical features, including central tendencies and
distribution properties, are extracted, providing insights into the dataset’s
shape and variability. Paralinguistic features are also analyzed to capture
non-verbal communication cues that are often atypical in children with ASD.
Acoustic features are examined to understand the fundamental frequency and
loudness, which can indicate differences in voice quality and pitch.
Linguistic features are assessed to capture vocabulary usage, sentence
structure, and language complexity, which can reveal differences in how
children with ASD use language compared to neurotypical children. The
resulting comprehensive set of features forms the basis for further analysis,
modeling, and a deep understanding of auditory characteristics associated with
ASD in children. We proceeded to employ a diverse set of classifiers of
machine and deep learning models. This strategic selection covers an in-depth
analysis of the subject and its application in the context of our problem
statement, allowing for a comprehensive evaluation of the extracted features
in the context of ASD detection. This critical phase of our research serves to
inform the subsequent sections, contributing towards enhancing early detection
strategies for ASD. While early language delays are frequently observed in
young children with ASD, they are not experienced by every individual.
Surprisingly, certain aspects of auditory processing appear to be heightened
in adolescents and adults diagnosed with ASD. This observation is supported by
an expanding body of research that explores the identification or
discrimination of pitch in both pure and complex tones. Additionally, studies
delve into the processing of pitch contour and rhythm in speech among
individuals with ASD. In addition to shedding light on the nature of
differences in non-word repetition performance, the imperative of matching
participants based on language ability becomes evident, given the pivotal role
that language level, specifically encompassing vocabulary size and associated
phonological knowledge, plays in determining proficiency in non-word
repetition. The ability to repeat non-words is influenced by their similarity
to real words and the phonotactic probability of non-words, suggesting that
stored lexical-phonological knowledge significantly influences the accuracy of
repetition. Various research groups posit that both vocabulary size and the
subsequent development of phonological knowledge are fundamental causal
factors influencing non-word repetition performance. Recognizing the critical
significance of language level and its related phonological knowledge in
shaping non-word repetition abilities, our study undertakes a comparative
analysis. We compare a group of kindergarten-age children diagnosed with ASD,
a majority of whom exhibit notable language delays, with a control group of
typically developing children carefully matched based on language ability.
Further, we explore the relationships between language proficiency or age, on
one hand, and non-word repetition accuracy, on the other, aspects that have
not been previously examined within the context of ASD. Our findings reveal a
substantial association between language level and non-word repetition
accuracy in ASD, mirroring established patterns in typical development, this
would provide additional compelling evidence for the presence of typical
phonological processes in this specific population.

##  2 Related Work

ASD is identified by repetitive behavior and challenges in social
communication and interaction across different situations Myers et al. (2007).
Linguistic deficits play a crucial role in both the diagnosis of ASD and early
identification, with initial referrals often based on concerns about language
delay Eigsti et al. (2011). Cohen et al. (2014). Also, social network analysis
showed that children with communication disorders exhibited more diverse
social interaction patterns than other children. Esposito et al. (2014)
distinguished ASD risk patterns in 15-month-old infants based on vocal
differences, including shorter cries and higher fundamental frequency. Young
children with ASD can be differentiated from typically developing peers based
on language content, such as the length of utterances and the number of nouns.
Research on the cross-context stability of social communication disturbances
in ASD is limited, partly due to a lack of efficient measurement methods. Luo
et al. (2016) used standard NLP methods to digitize and visualize descriptive
sentences, showing complex patterns distinguishing ASD. Studies in
neuroscience emphasize the significance of the corpus callosum and
intracranial brain volume for detecting ASD Sharif and Khan (2022).
Vocalization patterns across interactive and non-interactive contexts indicate
a strong predictive correlation for early communication development in ASD
Franklin et al. (2014). Spectrograms are used to model vocal patterns in the
speech of infants and children, revealing differences in energy distribution
Chittora and Patil (2016). Azizi (2015) demonstrated abnormal features in the
spectrographic analysis of children with ASD, assessing intonation patterns,
pitch, amplitude, duration, intensity, and tilt. The ASD population exhibits a
diverse range of acoustic speech patterns throughout language development, as
highlighted by studies conducted by Fusaroli et al. (2017); DePape et al.
(2012). In contrast to typically developing children, who embark on vocabulary
expansion around 24 months, children with ASD typically initiate this process
later Tager-Flusberg et al. (2005). The precise biological criteria and
origins of autism, along with its symptoms and fundamental characteristics,
remain unclear, as observed in studies by Herbert et al. (2005); Bonneh et al.
(2011) extensively documented supra segmental speech and extracted prosody
features in children with ASD, paving the way for subsequent studies.Various
dimensions of speech were explored in works by Quigley et al. (2016); Filipe
et al. (2014); Brisson et al. (2014); Sharda et al. (2010); Scharfstein et al.
(2011), disclosed insights into children and adolescents with ASD. Rahman et
al. (2022) provides a comprehensive review of automated ASD detection systems
based on human activity analysis (HAA), highlighting the effectiveness of
machine learning algorithms in identifying repetitive behaviors, abnormal
gait, and visual saliency patterns with high accuracy. Akhtar et al. (2024)
addresses the challenges of detecting ASD in code-switched conversations by
leveraging acoustic (A), linguistic (L), and paralinguistic (P) features,
achieving outstanding performance with a Transformer model through feature
fusion A+L then P. Ma et al. (2024) conducted a meta-analysis on the use of
natural speech prosody to distinguish individuals with (ASD from typically
developing (TD) individuals. The study found that pitch-related parameters
such as pitch mean, pitch range, pitch standard deviation, and pitch
variability were effective in identifying ASD, with significant effect sizes
and the ML models trained on these prosodic features.

##  3  CoSAm Corpus

In this study we gathered voice recordings from two groups: ASD and children
without ASD (the control group). Firstly, we conducted an comprehensive survey
with the purpose of identifying the institutions and caregivers who are
engaged in providing assistance and education to children with ASD. Once we
recognised these key stakeholders such as educational institutions, therapy
clinics, and government organizations that focus on supporting families
impacted by ASD, an NGO and its volunteers who had a comprehensive awareness
of the ASD landcape in the city, including a profound knowledge of the
different stages of autism among the children helped us to collect this
dataset.

###  3.1 Selecting Participant

To be part of the group with ASD, a child had to have a clinicians diagnosis
of ASD and a part of autistic welfare center. On the other hand, the healthy
group, it was confirmed from the clinicians, that they do not exhibit any
symptoms and are a part of conventional family settings. Additionally, extra
precautions were taken to confirm that both the children with ASD and the non-
ASD had no neurological disorders or any record of head injuries. Notably, all
the selected children diagnosed with ASD consistently exhibited unconventional
speech patterns throughout our data collection phase.

###  3.2 Recording Audios

We managed to collect audios from a small number of samples or population for
our statistical and acoustics analysis. The ASD-affected and the healthy
children were instructed to follow the provided question script (The scripts
can be provided on request). The script consists of 15 questions in both
English and Hindi, and the students were trained to follow and respond to the
script. For preparing our script, we used visual cues from the primary
learning textbooks containing visual depictions illustrating vegetables,
animals, and other subject matter aligned with their curriculum. The reason
being children’s comfort and familiarity with the questions asked in the
session. The speech samples were captured in stereo mode, employing four
channels. The audio recordings of the speech sound samples were saved in .wav
format with a sampling rate of 48 kHz. Speech samples were recorded in an
empty room, ensuring the absence of background noise of any objects that could
generate noise or distract the children during the recording sessions.
Moreover, to minimize background noise, all fans were turned off throughout
the recording sessions. The recording equipment was strategically positioned
to prevent any echoes caused by the children’s voices during the sessions. All
children were instructed to sit in a chair for the recording of speech
samples, and the recording equipment was set up on a table for optimal
placement.

###  3.3  CoSAm: Statistics

The dataset consists of 61 voice recordings from 30 children diagnosed with
ASD from mild to moderate condition and 31 from neuro-typical children, aged
from 3 and 13 years. All the voice recordings are in code-mixing scenario of
English and Hindi. CoSAm corpus has total 159.75 minutes of voice recordings
with total segmented sample of 151 and 62 for ASD and healthy control
respectively. We would release our corpus based on the request and strictly
for research purpose. Additionally, apart from the voice recordings we will
also release the corresponding transcripts from Whisper Radford et al. (2022)
to further advance the future research in this domain. We have provided the
detailed statistics of our corpus in Table 1.

All the recordings were obtained in a regulated setting to reduce ambient
noise and guarantee uniformity in audio quality. The recordings were captured
in a controlled environment to minimize background noise and ensure
consistency in audio quality. Each recording session was designed to encourage
natural vocal expressions, with prompts and activities tailored to elicit a
range of vocalizations and speech patterns. We present visualizations of the
gender distribution, class/category distribution, age distribution, and the
trends in symptom proportions across age and gender within CoSAm. Figure:1(a)
shows the gender distribution, indicating that 62.3% are male (M) and 37.7%
are female (F) audios. Figure: 1(b) shows gender distribution for healthy
control or non-autism and autism and are represented with blue and red bar
respectively. It is evident from this figure that number of males are higher
compared to the number of females across the classes. All the audios, are
recorded for children from 6-13 years and age distribution box plot is shown
in Figure:1(c). It shows the median, the interquartile range (the middle 50%
of the data), and the minimum and maximum values, giving insights into the
spread and skewness of the age distribution. Finally in Figure: 1(d), we
represent the proportion of symptoms across different age groups, with the
blue line representing males and the red line representing females. The plot
shows how the proportion of symptoms varies as a function of age for both
genders.

Table 1: Dataset Details Attribute | Group | Male | Female | Total
---|---|---|---|---
Number of Children | ASD | 18 | 12 | 30
Non-ASD | 20 | 11 | 31
Age (years) | ASD | 06 to 13
Non-ASD | 06 to 13
Data Duration (min) | ASD | 105.76
Non-ASD | 53.99
Segmented Samples | ASD | 151
Non-ASD | 62
Native Language | ASD | Codemixed
Non-ASD | Codemixed
Level | ASD | Mild/Moderate
Non-ASD | Normal

###  3.4 Handling Ethical Issues and Challenges

For our corpus, we had a prior consent of all the participants, caregivers or
legal guardians of the children who participated in our study. We further de-
identify our dataset, and ensure full anonymity and well-being of our
participants. Later, we followed an inclusion and exclusion criteria, by
involving an healthcare professional for labelling our CoSAm corpus. Curating
CoSAm was challenging task due to the atypical or unpredicted behavior of the
children with ASD. Throughout the recording sessions, a few rather hostile
kids sometimes began yelling. Most of the time, they made an effort to speak
the words loudly. Also, several kids were very nervous during the recording
sessions. It was quite difficult to keep the kids’ emotions in check
throughout the recording sessions and capture the speech samples with perfect
clarity in such circumstances. Therefore, we took several intervals between
recording sessions. We also make engaged children in a few enjoyable
activities. Another significant challenge was that many selected children with
ASD would, at times, utter extra words in between recording sessions, adding
complexity to the recording process. To address this issue, we used the
PRAAT111https://praat.en.softonic.com/ tool to manually remove the extra words
and gaps in between words, thereby ensuring the accuracy of the recorded
speech samples. Our collected corpus is novel and to the best of our knowledge
first autistic children corpus with code-mixing scenarios. Primarily, both of
our databases feature exclusively native English and Hindi speakers,
encompassing all selected speakers, including children with ASD and non-ASD
children.

###  3.5 Statistical Analysis

Table 2: Statistical analysis of audio recordings | ENTIRE DATASET | TRAIN DATASET
---|---|---
FEATURE | NORMAL | AUTISM | NORMAL | AUTISM
FEMALE
Meanf0 | 35.42±3.23 | 35.31±2.39 | 35.27±3.13 | 35.41±2.48
StdDevf0 | 0.17±0.06 | 0.13±0.04 | 0.17±0.06 | 0.14±0.05
Mean Loudness | 0.95±0.77 | 0.50±0.20 | 0.88±0.73 | 0.49±0.25
Stddev Loudness | 0.77±0.13 | 0.91±0.13 | 0.76±0.13 | 0.92±0.13
Localshimmer Mean | 1.15±0.13 | 1.06±0.18 | 1.71±0.12 | 1.04±0.18
LocaldBshimmer | 0.80±0.11 | 0.89±0.14 | 0.78±0.09 | 0.89±0.15
HNRdB | 6.83±1.73 | 8.87±1.51 | 6.77±1.62 | 8.97±1.48
F1-Frequency | 602.13±79.21 | 686.95±41.67 | 597.61±81.85 | 685.33±43.85
Local Jitter Mean | 0.04±0.01 | 0.02±0.006 | 0.04±0.01 | 0.02±0.005
Local Jitter stddev | 1.89±0.28 | 1.79±0.36 | 1.88±0.27 | 1.81±0.38
LongRelf0mean | 8.70±2.71 | 8.35±2.88 | 8.53±2.59 | 8.22±2.72
LongRelf0stddev | 1.71±1.004 | 1.46±0.588 | 1.64±0.69 | 1.48±0.58
MALE
Meanf0 | 33.13±3.28 | 35.30±2.10 | 32.68±3.41 | 35.33±2.01
StdDevf0 | 0.21±0.07 | 0.15±0.04 | 0.23±0.07 | 0.15±0.04
Mean Loudness | 0.62±0.29 | 0.53±0.17 | 0.65±0.32 | 0.53±0.17
Stddev Loudness | 0.72±0.11 | 0.86±0.13 | 0.72±0.13 | 0.87±0.12
Localshimmer Mean | 1.16±0.13 | 1.06±0.14 | 1.81±0.13 | 1.06±0.14
LocaldBshimmer | 0.85±0.08 | 0.89±0.11 | 0.86±0.08 | 0.89±0.12
HNRdB | 6.49±2.14 | 9.01±1.34 | 6.26±2.32 | 8.98±1.39
F1-Frequency | 577.81±44.03 | 687.62±39.29 | 585.11±43.49 | 685.93±38.30
Local Jitter Mean | 0.04±0.01 | 0.02±0.005 | 0.04±0.02 | 0.02±0.004
Local Jitter stddev | 1.77±0.33 | 1.70±0.30 | 1.71±0.34 | 1.71±0.30
Longrelf0mean | 8.42±2.93 | 2.25±0.65 | 7.96±2.99 | 8.51±2.844
LongRelf0stddev | 2.05±1.81 | 0.29±0.19 | 2.27±1.97 | 1.47±0.55

Table 2 presents analysis of audio recordings and it reveals distinct
differences between individuals classified as Normal and those with Autism,
across both the entire data set and the train data set, further classified by
gender. For females, the mean fundamental frequency (Meanf0) is nearly
identical between the normal (35.42±3.23 Hz) and autism (35.31±2.39 Hz) groups
in the entire dataset, and similar trends are observed in the train dataset
(normal: 35.27±3.13 Hz, autism: 35.41±2.48 Hz). The variability (StdDevf0) is
slightly lower in the autism group (0.13±0.04 Hz) compared to the normal group
(0.17±0.06 Hz) in the entire dataset. Mean loudness is significantly lower in
the autism group (0.50±0.20) compared to the normal group (0.95±0.77) in the
entire dataset, and this pattern remains in the train dataset (autism:
0.49±0.25, normal: 0.88±0.73). Variability in loudness is higher in the autism
group (0.91±0.13) than in the normal group (0.77±0.13). Females with autism
also have a higher harmonics-to-noise ratio (HNRdB) (8.87±1.51) compared to
the normal group (6.83±1.73). On the other hand, for males, the autism group
has a higher mean fundamental frequency (35.30±2.10 Hz) than the normal group
(33.13±3.28 Hz) in the entire dataset, with similar results in the train
dataset (autism: 35.33±2.01 Hz, normal: 32.68±3.41 Hz). The variability in
fundamental frequency (StdDevf0) is lower in the autism group (0.15±0.04 Hz)
compared to the normal group (0.21±0.07 Hz). Mean loudness is comparable
between the groups (autism: 0.53±0.17, normal: 0.62±0.29), but the autism
group shows higher variability in loudness (0.86±0.13) than the normal group
(0.72±0.11). The autism group also demonstrates a higher HNRdB (9.01±1.34)
compared to the normal group (6.49±2.14). These findings suggest that
individuals with autism exhibit unique acoustic features, including lower mean
loudness, higher variability in loudness, and a higher harmonics-to-noise
ratio, which could be instrumental for diagnostic or therapeutic applications.

(a)

01$0$$5$$10$$15$$20$ClassCountMalesFemales (b)

01$7$$8$$9$$10$$11$$12$$13$ClassAge (c)

$7$$8$$9$$10$$11$$12$$13$$14$$1.3$$1.4$$1.5$$1.6$$1.7$$1.8$$1.9$$2$AgeProportion
of AbnormalFM (d)

Figure 1: Distribution of dataset on autistic children. (a) shows the gender
distribution, with 62.9% male and 37.1% female. (b) Displays the
class/category distribution without sgemented samples, with higher
representation of males compared to females across the different classes. (c)
Box plot illustrating the age distribution, providing insights into the
median, interquartile range, and minimum and maximum values. (d) Line plot
depicting the proportion of symptoms across different age groups for males and
females; 0 and 1 in Figure 1(b) and 1(c) represents healthy control and
autistic children; F and M stand for Female and Male, respectively

###  3.6 Feature Extraction

Our study emphasizes the comprehensive analysis of speech samples through a
multi-faceted feature extraction approach that includes Acoustic, Linguistic,
and Paralinguistic features using the XLSR-53, mBERT, and TRILLsson models
respectively. Here, we describe each feature type and the integration methods
used in our experiments.

Acoustic features: We leverage XLSR-53 Grosman (2021) for taking out the
acoustic features. It analyzes the raw audio signals to extract strong
features that capture speech’s time and frequency aspects, resulting in
comprehensive representations. The decision to select XLSR-53 was based on its
demonstrated effectiveness in various applications of multilingual speech
processing, attributed to its pre-training in diverse data from 53 different
languages. This makes it especially appropriate for studying the speech of
children with diverse linguistic origins.

Linguistic Features: In order to analyze the linguistic aspects of speech, we
used two approaches for extracting the required features. Firstly, we use
prominent handcrafted features primarily included in code-mixed language
processing research using NLTK library, such as average word length, average
sentence length, speech rate, language distribution (English/Hindi), and code-
switching patterns. Secondly, we considered multilingual BERT (mBERT) Devlin
et al. (2018). It is pre-trained on 104 languages using masked language
modeling and next sentence prediction objective. This approach enabled us to
examine the subtle differences in meaning and sentence structure seen in the
language of children with and without ASD. Joint training on multiple
languages provides contextual embeddings that captures lexical diversity,
semantic coherence, and syntactic complexity of several different languages.

Paralinguistic Features: We use TRILLsson Shor and Venugopalan (2022) for
extracting paralinguistic features. TRILLsson is distilled from the larger
CAP12 model. It retains high performance while significantly reducing the
model size to 1%-15% of CAP12, which has 600 million parameters and is trained
on over 900 million hours of YouTube audio. TRILLsson learns the non-verbal
components of speech such as intonation, stress, and rhythm. These features
are particularly salient in differentiating the speech patterns of children
with ASD, who often exhibit atypical paralinguistic characteristics.

##  4 Modeling

We divide our modeling section into three subsections: individual
representation modeling, fusion with concatenation modeling, and hierarchical
fusion modeling. The performance metrics used for evaluation included accuracy
and macro F1 score. We employed a 5-fold cross-validation approach, where 4
folds were used for training and 1 fold for testing.

Individual representation modeling: We build models for linguistic (L),
acoustic (A), and paralinguistic (P) features separately. Figure 2(a)
illustrates the architecture for individual modeling. We experiment with both
classical ML algorithms such as SVM, RF, kNN, NB, and DT and deep learning
techniques such as RNN, CNN, and Transformer. For classical ML algorithms, we
trained the models with the default hyperparameters as set in Scikit-learn
library. For DL models, for the RNN approach, we employed a single-layer RNN
architecture each consisting of 50 hidden cells. Next for CNN, each feature
type was processed through a 1D-CNN layer with 64 filters of size 3, followed
by a max-pooling layer with a pool size of 2. For the Transformer approach, we
utilized the same initial processing as the CNN, incorporating a 1D-CNN layer
with 64 filters of size 3 and a max-pooling layer, followed by a transformer
encoder with number of attention heads as 4. We added dense layers with 128
and 64 neurons respectively as the classifier with the output layer
representing the class ASD and non-ASD.

Fusion with Concatenation: We combine the individual representations through
concatenation modeling to enhance the model’s performance. Various fusion
scenarios were explored: L + A, A + P, L + P, and A + L + P. For each
combination, we evaluated both ML and DL models. Figure 2(b) depicts the
concatenation modeling with the transformer-based approach. The modeling here
follows the same architectural details as in the above individual
representation modeling. For classical ML algorithms, the parameters are used
as default. We use concatenation as the fusion mechanism, followed by the
dense layers of 128 and 64 neurons as the classifier.

Hierarchical Feature Fusion: The final phase involved a hierarchical feature
fusion approach to further improve the model’s performance. We present our
novel hierarchical feature fusion method in Figures 2(c),2(d),2(e). Here, we
also keep the modeling details same as the concatenation modeling approach
above. However, we use a dense layer of 120 neurons in the hierarchical fusion
part and followed by using concatenation for fusion. The hierarchical fusion
approach aimed to optimize the combined performance of features, leveraging
the strengths of each modality to enhance the model’s capabilities.

For all deep learning models, including RNN, CNN, and Transformer, we used
Cross-entropy as loss function, adam as optimizer with a learning rate of
0.001, batch size of 32, and trained each model for 100 epochs. Early stopping
and dropout were used to reduce overfitting.

(a)

(b)

(c)

(d)

(e)

Figure 2: Architecture: (a) Individual Modeling; (b) Feature Fusion Through
Concatenation; (c,d,e) Hierarchical Feature Fusion; Here, the architectural
diagrams are only given for the transformer modeling approach

##  5 Experiments

In this section, we brief our experiments, for our proposed hierarchical
feature fusion technique in identifying ASD by analyzing children’s speech. In
our experiments, we examine the idea that speech patterns might be used as
dependable biomarkers for the early identification of autism.

Figure 3: Confusion Matrices: (a) Linguistic (Transformer), (b) Acoustic
(Transformer), (c) Linguistic-mBERT (Transformer), (d) A+L then P
(Transformer)

|  | Predicted
---|---|---
Actual |  | 0 | 1
0 | 12 | 1
1 | 3 | 27
(a)

|  | Predicted
---|---|---
Actual |  | 0 | 1
0 | 11 | 2
1 | 1 | 29
(b)

|  | Predicted
---|---|---
Actual |  | 0 | 1
0 | 11 | 2
1 | 0 | 30
(c)

|  | Predicted
---|---|---
Actual |  | 0 | 1
0 | 13 | 0
1 | 1 | 29
(d)

###  5.1 Experiment Setup

All of our research experiments were conducted using Python 3.9. We use a
range of specialized libraries, such as Librosa for audio processing, Hugging
Face’s Transformers for accessing pre-trained models, TensorFlow and Keras for
developing deep learning models, and SciPy, NumPy, and Pandas for analyzing
and manipulating data. For our experimental setup, we used a high-performance
computing setup featuring a Ryzen 7 processor, NVIDIA GeForce RTX 3050 GPU, 32
GB RAM, and 1 TB SSD to handle our extensive data processing and model
training needs. This configuration guarantees the ability to replicate our
findings and makes it easier to make additional progress in the early
identification of ASD using speech analysis by using the strong computational
framework described.

###  5.2 Experimental Results

For the evaluation of our proposed hierarchical feature fusion method in
detecting autism using code-switched speech data, we conducted a series of
experiments using various machine learning and deep learning models.
As a baseline model, we used MFCC features. The results are presented in Table
3, highlighting the accuracy and macro F1 score. SVM and RF demonstrated
moderate performance, achieving accuracy around 81-82% and macro F1 scores
around 68-69%. In contrast, KNN and NB exhibited lower accuracy of
approximately 58-59% and macro F1 scores of around 57-58%, indicating
challenges in effectively modeling the data. DT performed better than KNN and
NB, but still fell short of SVM and RF, with an accuracy of around 70% and a
macro F1 score of approximately 62%. The RNN achieved relatively high accuracy
at 71.16%, but its macro F1 score was notably lower at 41.51%, suggesting
issues with precision and recall. The best performer was the transformer
model, achieving the highest accuracy of 89.30% and a strong macro F1 score of
86.95%, indicating its effectiveness in capturing complex patterns in the
data. Below we give feature-wise modeling results.

Individual Feature Modeling: We evaluated the performance of different models
utilizing individual representations of A, L and P features. These results are
shown in Table 4. SVM showed the best performance among traditional ML models,
achieving an accuracy of 86.98% and a macro F1 score of 81.44%. Similarly, RF
achieved an accuracy of 85.12% and a macro F1 score of 80.80%. KNN showed
strong performance, especially with A and P features, with an accuracy of
81.40% and a macro F1 score of 75.27%. NB and DT exhibited decent but
comparatively lower accuracy, achieving 80.47% and 79.07%, respectively, with
macro F1 scores of 74.19% and 75.38%. The transformer-based architectures
outperformed other models, particularly with L and A features, achieving an
accuracy of 87.90%, with strong macro F1 scores above 83%.
Fusion of Feature Representations Through Concatenation: The results are
displayed in Table 5. We considered four fusion scenarios: L + A, A + P, L +
P, and A + L + P. For each model, we reported accuracy and macro F1
performance. SVM consistently achieved high accuracy and macro F1 scores in
most fusion scenarios among classical ML algorithms. Among the DL techniques,
transformer models achieved peak performance of 98.13% accuracy and 98.00%
macro F1 score in the L + P fusion. RF performed competitively but fell
slightly behind SVM in most scenarios, reaching a maximum of 96.88% accuracy
and 96.10% macro F1 score in the A + L + P fusion. KNN demonstrated strong
performance across all fusion scenarios, particularly in the A + L + P fusion
with 97.50% accuracy and 96.93% macro F1 score. NB and DT showed decent
performance but generally fell short compared to SVM, RF, and KNN. RNN and CNN
achieved moderate performance across fusion scenarios but showed lower
performance compared to SVM, RF, KNN, and Transformer-based models. The
Transformer consistently outperformed other models in most fusion scenarios,
demonstrating particularly high accuracy and macro F1 scores in the A + P and
L + P combinations. The Transformer achieved a peak performance of 98.31%
accuracy and 98.00% macro F1 score in the A + L + P fusion where mBERT was
used for linguistic features, while the Transformer reached 98.14% accuracy
and 98.00% macro F1 score in the A + P fusion and 98.13% accuracy and 98.00%
macro F1 score in the L + P fusion.
Hierarchical Feature Fusion: Here, we fused these features in hierarchical
order, exploring how different modalities impact performance when combined in
specific sequences. Table 8 reports the performance of the transformer model
based on accuracy and macro F1 score. Firstly, the A+L THEN P combination
yielded impressive results, achieving 98.75% accuracy and a macro F1 score of
97.19%. This suggests that incorporating P and A features before L features
leads to exceptional performance. Conversely, the A+L THEN P combination
achieved slightly lower accuracy at 95.81% and a macro F1 score of 94.00%.
While still performing well, this sequence fell slightly short compared to the
previous combination. The L+P THEN A combination emerged as the best
performing, with an accuracy of 98.60% and a macro F1 score of 97.00%. In a
hierarchical context specifically focusing on L features processed by the
mBERT model, the P+A THEN L combination achieved an accuracy of 97.50% and a
macro F1 score of 97.00%. The A+L THEN P combination stood out with the
highest accuracy of 98.75% while maintaining a macro F1 score of 97.00%. The
L+P THEN A combination achieved an accuracy of 96.56% and a macro F1 score of
96.00%, indicating slightly lower performance compared to the other
combinations. These findings underscore the significance of the sequence in
which modalities are combined, with certain combinations yielding notably
higher accuracy and F1 scores. The results highlight the potential for
optimizing model performance through strategic hierarchical feature
integration, particularly in tasks related to linguistic analysis and
representation modeling. The transformer consistently outperformed, achieving
the highest accuracy and macro F1 scores across feature combinations. The
performance of DL models, particularly the transformer, was superior in
capturing complex patterns present in different feature combinations.

Table 3: MFCC Performance Analysis; Scores presented are average of 5 folds; Scores are in % Model | Accuracy | F1
---|---|---
SVM | 81.86 | 69.34
RF | 81.79 | 68.13
KNN | 58.14 | 57.65
NB | 59.07 | 57.89
DT | 70.23 | 61.79
RNN | 71.16 | 41.51
CNN | 75.81 | 67.95
TRANSFORMER | 89.30 | 86.95
Table 4: Performance Analysis: Individual Representation; Scores presented are average of 5 folds; Scores are in % Model | LINGUISTIC | ACOUSTIC | PARALINGUISTIC
---|---|---|---
| Accuracy | F1 | Accuracy | F1 | Accuracy | F1
SVM | 86.98 | 81.44 | 90.23 | 86.66 | 97.38 | 96.84
RF | 85.12 | 80.80 | 87.91 | 83.30 | 94.77 | 93.50
KNN | 81.40 | 75.27 | 89.30 | 87.72 | 97.57 | 97.11
NB | 80.47 | 74.19 | 84.65 | 81.27 | 92.34 | 90.41
DT | 79.07 | 75.38 | 76.28 | 72.71 | 89.72 | 87.18
RNN | 85.12 | 78.90 | 86.51 | 81.42 | 83.72 | 76.82
CNN | 87.44 | 80.49 | 89.30 | 86.34 | 98.13 | 97.37
TRANSFORMER | 87.90 | 83.50 | 93.48 | 91.50 | 97.67 | 96.50
Table 5: Performance Analysis: Fusion with Concatenation; Scores presented are average of 5 folds; For SVM,RF,KNN,NB,DT Fusion of the features is done before passing through the models; Scores are in % Model | L+A | A+P | L+P | A+L+P
---|---|---|---|---
| Accuracy | F1 | Accuracy | F1 | Accuracy | F1 | Accuracy | F1
SVM | 90.23 | 86.66 | 97.08 | 96.36 | 97.50 | 96.91 | 96.12 | 94.66
RF | 86.98 | 82.55 | 86.88 | 86.11 | 95.63 | 94.52 | 96.88 | 96.10
KNN | 89.77 | 88.17 | 97.50 | 96.89 | 96.09 | 95.24 | 97.50 | 96.93
NB | 84.65 | 80.81 | 92.08 | 86.98 | 94.84 | 93.35 | 92.81 | 90.55
DT | 77.21 | 73.39 | 91.23 | 88.79 | 89.53 | 87.06 | 89.38 | 86.51
RNN | 83.72 | 77.17 | 82.33 | 74.40 | 84.19 | 77.87 | 90.23 | 86.27
CNN | 91.16 | 88.83 | 97.67 | 97.37 | 96.27 | 95.16 | 97.20 | 96.86
TRANSFORMER | 92.09 | 90.50 | 98.14 | 98.00 | 98.13 | 98.00 | 96.74 | 97.00
Table 6: Performance Analysis: Hierarchical; Scores presented are average of 5 folds; The model used here is transformer; Scores are in % Model/Combination |  Accuracy |  F1
---|---|---
P+A THEN L |  98.13 |  96.00
A+L THEN P |  95.81 |  94.00
L+P THEN A |  98.60 |  97.00
Table 7: Combined Performance Analysis; Scores presented are average of 5 folds; mBERT is used for linguistic (L) feature extraction; Scores are in % Model | L | L+A | L+P | A+L+P
---|---|---|---|---
| Accuracy | F1 | Accuracy | F1 | Accuracy | F1 | Accuracy | F1
SVM | 89.69 | 85.28 | 91.41 | 87.91 | 96.64 | 95.81 | 95.89 | 94.75
RF | 87.34 | 82.17 | 88.91 | 84.55 | 94.39 | 93.02 | 92.90 | 90.83
KNN | 86.56 | 80.53 | 90.00 | 86.10 | 96.82 | 96.18 | 96.45 | 95.72
NB | 86.25 | 81.89 | 89.06 | 85.21 | 94.58 | 92.99 | 93.27 | 91.37
DT | 85.00 | 81.32 | 85.78 | 81.89 | 89.91 | 87.69 | 89.53 | 87.07
RNN | 81.82 | 76.34 | 78.14 | 68.54 | 82.73 | 78.85 | 88.18 | 84.34
CNN | 94.41 | 93.03 | 95.34 | 94.14 | 96.74 | 96.00 | 97.20 | 96.81
TRANSFORMER | 95.81 | 95.00 | 94.01 | 92.50 | 98.14 | 98.00 | 98.31 | 98.00
Table 8: Hierarchical Performance Analysis; mBERT is used for linguistic (L) feature extraction; The model used here is transformer; Scores are in % Model | Accuracy | Macro F1
---|---|---
P+A THEN L | 97.50 | 97.10
A+L THEN P | 98.75 | 97.19
L+P THEN A | 96.56 | 96.03

##  6 Discussion

We present a novel hierarchical feature fusion method for early detection of
Autism Spectrum Disorder (ASD) using the CoSAm code-switched speech corpus
(English and Hindi). By integrating A, L, and P features through transformer
encoders, we significantly enhanced the robustness and accuracy of ASD
identification. The hierarchical fusion technique demonstrated superior
performance compared to existing methods, achieving an accuracy of 98.75% when
combining A and L features, followed by P features. The confusion matrices in
figure 3 provide a detailed breakdown of our model’s performance. Figure 3(d)
presents the confusion matrix for the hierarchical fusion of A+L and then P
using a transformer, which shows the best performance. This model achieved a
high true positive rate and minimal false positives and false negatives,
underscoring its effectiveness.

This finding underscores the critical importance of the feature integration
sequence, highlighting the role of modality order in optimizing model
performance. The success of our hierarchical approach is particularly notable
when compared to the performance of models trained on MFCC features which were
preferred features by previous prolific works on autism detection Hu et al.
(2024), Kim et al. (2021), Chowdhury et al. (2024). Our experiments showcases
that effective selection of features and the method of fusing the features
leads to better ASD detection.

##  7 Conclusion

We introduce CoSAm, a unique dataset featuring code-switched speech from
bilingual children with ASD. The study introduces a novel hierarchical feature
fusion method designed to enhance the early detection of ASD in children
through the analysis of code-switched speech (English and Hindi). By
leveraging advanced speech processing techniques, our approach effectively
integrates A, L, and P features, resulting in a robust and accurate
classification model. Our experimental results demonstrate that the sequence
in which features are fused significantly impacts model performance.
Specifically, we evaluated several hierarchical fusion sequences: A and L
features followed by P features (A+L then P), P and A features followed by L
features (P+A then L), and L and P features followed by A features (L+P then
A). Among these combinations, the A+L then P fusion achieved the topmost
accuracy of 98.75% and a macro F1 score of 97.19%, highlighting the
effectiveness of this sequence in capturing complex speech patterns associated
with ASD. On the other hand, the P+A then L combination achieved slightly
lower performance, with an accuracy of 97.50% and a macro F1 score of 97.10%.
The L+P then A combination, while still effective, demonstrated an accuracy of
96.56% and a macro F1 score of 96.03%. These findings underscore the critical
importance of modality order in feature integration, showcasing the
effectiveness of our proposed method in optimizing model performance.
Our results highlight the potential of sophisticated signal processing
techniques in the early detection of neurological conditions. The hierarchical
fusion technique demonstrated superior performance compared to existing
methods, achieving outstanding accuracy when combining A and L features,
followed by P features. This underscores the significance of the sequence in
which features are integrated and provides valuable insights into the speech
patterns of bilingual children with ASD.
Future Work: We will further explore refinements in feature integration
strategies and extend this approach to other datasets and languages to
validate and enhance its generalizability.

##  8 Source Code and Data Availability Statement

The dataset generated and analyzed during research is available from the
corresponding author on reasonable request for research purposes. The source
code used in this study will be available at: https://github.com/gir-
ish/Autism

## References

  * McCann and Peppé (2003) J. McCann, S. Peppé,  Prosody in autism spectrum disorders: a critical review,  International Journal of Language & Communication Disorders 38 (2003) 325–350.
  * Diehl et al. (2009) J. J. Diehl, D. Watson, L. Bennetto, J. McDonough, C. Gunlogson,  An acoustic analysis of prosody in high-functioning autism,  Applied Psycholinguistics 30 (2009) 385–404.
  * Kjelgaard and Tager-Flusberg (2001) M. M. Kjelgaard, H. Tager-Flusberg,  An investigation of language impairment in autism: Implications for genetic subgroups,  Language and cognitive processes 16 (2001) 287–308.
  * Myers et al. (2007) S. M. Myers, C. P. Johnson, C. on Children with Disabilities,  Management of children with autism spectrum disorders,  Pediatrics 120 (2007) 1162–1182.
  * Eigsti et al. (2011) I.-M. Eigsti, A. B. de Marchena, J. M. Schuh, E. Kelley,  Language acquisition in autism spectrum disorders: A developmental review,  Research in Autism Spectrum Disorders 5 (2011) 681–691.
  * Cohen et al. (2014) I. L. Cohen, J. M. Gardner, B. Z. Karmel, S.-Y. Kim,  Rating scale measures are associated with noldus ethovision-xt video tracking of behaviors of children on the autism spectrum,  Molecular autism 5 (2014) 1–17.
  * Esposito et al. (2014) G. Esposito, M. del Carmen Rostagno, P. Venuti, J. D. Haltigan, D. S. Messinger,  Brief report: Atypical expression of distress during the separation phase of the strange situation procedure in infant siblings at high risk for asd,  Journal of Autism and Developmental Disorders 44 (2014) 975–980.
  * Luo et al. (2016) S. X. Luo, J. A. Shinall, B. S. Peterson, A. J. Gerber,  Semantic mapping reveals distinct patterns in descriptions of social relations in adults with autism spectrum disorder,  Autism Research 9 (2016) 846–853.
  * Sharif and Khan (2022) H. Sharif, R. A. Khan,  A novel machine learning based framework for detection of autism spectrum disorder (asd),  Applied Artificial Intelligence 36 (2022) 2004655.
  * Franklin et al. (2014) B. Franklin, A. S. Warlaumont, D. Messinger, E. Bene, S. Nathani Iyer, C.-C. Lee, B. Lambert, D. K. Oller,  Effects of parental interaction on infant vocalization rate, variability and vocal type,  Language Learning and Development 10 (2014) 279–296.
  * Chittora and Patil (2016) A. Chittora, H. A. Patil,  Spectral analysis of infant cries and adult speech,  International Journal of Speech Technology 19 (2016) 841–856.
  * Azizi (2015) Z. Azizi,  The acoustic survey of intonation in autism spectrum disorder,  The Journal of the Acoustical Society of America 137 (2015) 2207–2207.
  * Fusaroli et al. (2017) R. Fusaroli, A. Lambrechts, D. Bang, D. M. Bowler, S. B. Gaigg,  Is voice a marker for autism spectrum disorder? a systematic review and meta-analysis,  Autism Research 10 (2017) 384–407.
  * DePape et al. (2012) A.-M. R. DePape, A. Chen, G. B. Hall, L. J. Trainor,  Use of prosody and information structure in high functioning adults with autism in relation to language ability,  Frontiers in psychology 3 (2012) 21096.
  * Tager-Flusberg et al. (2005) H. Tager-Flusberg, R. Paul, C. Lord,  Language and communication in autism,  Handbook of autism and pervasive developmental disorders 1 (2005) 335–364.
  * Herbert et al. (2005) M. R. Herbert, et al.,  Autism: a brain disorder or a disorder that affects the brain,  Clinical Neuropsychiatry 2 (2005) 354–379.
  * Bonneh et al. (2011) Y. S. Bonneh, Y. Levanon, O. Dean-Pardo, L. Lossos, Y. Adini,  Abnormal speech spectrum and increased pitch variability in young autistic children,  Frontiers in human neuroscience 4 (2011) 237.
  * Quigley et al. (2016) J. Quigley, S. McNally, S. Lawson,  Prosodic patterns in interaction of low-risk and at-risk-of-autism spectrum disorders infants and their mothers at 12 and 18 months,  Language Learning and Development 12 (2016) 295–310.
  * Filipe et al. (2014) M. G. Filipe, S. Frota, S. L. Castro, S. G. Vicente,  Atypical prosody in asperger syndrome: Perceptual and acoustic measurements,  Journal of autism and developmental disorders 44 (2014) 1972–1981.
  * Brisson et al. (2014) J. Brisson, K. Martel, J. Serres, S. Sirois, J.-L. Adrien,  Acoustic analysis of oral productions of infants later diagnosed with autism and their mother,  Infant mental health journal 35 (2014) 285–295.
  * Sharda et al. (2010) M. Sharda, T. P. Subhadra, S. Sahay, C. Nagaraja, L. Singh, R. Mishra, A. Sen, N. Singhal, D. Erickson, N. C. Singh,  Sounds of melody—pitch patterns of speech in autism,  Neuroscience letters 478 (2010) 42–45.
  * Scharfstein et al. (2011) L. A. Scharfstein, D. C. Beidel, V. K. Sims, L. Rendon Finnell,  Social skills deficits and vocal characteristics of children with social phobia or asperger’s disorder: A comparative study,  Journal of abnormal child psychology 39 (2011) 865–875.
  * Rahman et al. (2022) S. Rahman, S. F. Ahmed, O. Shahid, M. A. Arrafi, M. A. R. Ahad,  Automated detection approaches to autism spectrum disorder based on human activity analysis: A review,  Cognitive Computation 14 (2022) 1773–1800.
  * Akhtar et al. (2024) M. M. Akhtar, O. C. Phukan, M. Singh, et al.,  Neuro: An application for code-switched autism detection in children,  arXiv preprint arXiv:2406.03514 (2024).
  * Ma et al. (2024) W. Ma, L. Xu, H. Zhang, S. Zhang,  Can natural speech prosody distinguish autism spectrum disorders? a meta-analysis,  Behavioral Sciences 14 (2024) 90.
  * Radford et al. (2022) A. Radford, J. W. Kim, T. Xu, G. Brockman, C. McLeavey, I. Sutskever, Robust speech recognition via large-scale weak supervision, 2022. URL: https://arxiv.org/abs/2212.04356. doi:10.48550/ARXIV.2212.04356.
  * Grosman (2021) J. Grosman, Fine-tuned XLSR-53 large model for speech recognition in English, https://huggingface.co/jonatasgrosman/wav2vec2-large-xlsr-53-english, 2021.
  * Devlin et al. (2018) J. Devlin, M. Chang, K. Lee, K. Toutanova,  BERT: pre-training of deep bidirectional transformers for language understanding,  CoRR abs/1810.04805 (2018).
  * Shor and Venugopalan (2022) J. Shor, S. Venugopalan,  Trillsson: Distilled universal paralinguistic speech representations,  arXiv preprint arXiv:2203.00236 (2022).
  * Hu et al. (2024) C. Hu, J. Thrasher, W. Li, M. Ruan, X. Yu, L. K. Paul, S. Wang, X. Li,  Exploring speech pattern disorders in autism using machine learning,  arXiv preprint arXiv:2405.05126 (2024).
  * Kim et al. (2021) Y.-K. Kim, R. Lahiri, M. Nasir, S. H. Kim, S. Bishop, C. Lord, S. S. Narayanan,  Analyzing short term dynamic speech features for understanding behavioral traits of children with autism spectrum disorder.,  in: Interspeech, 2021, pp. 2916–2920.
  * Chowdhury et al. (2024) T. Chowdhury, V. Romero, A. Stent,  Parameter selection for analyzing conversations with autism spectrum disorder,  arXiv preprint arXiv:2401.09717 (2024).


# Band engineered bilayer Haldane model: Evidence of multiple topological
phase transitions

Sayan Mondal  Department of Physics, Indian Institute of Technology Guwahati,
Guwahati 781039, Assam, India  Saurabh Basu  Department of Physics, Indian
Institute of Technology Guwahati, Guwahati 781039, Assam, India

###### Abstract

We have studied the evolution of the topological properties of a band-
engineered AB-stacked bilayer honeycomb structure in the presence of a Haldane
flux. Without a Haldane flux, band engineering makes the band touching points
(the so called Dirac points) move towards each other and eventually merge into
one at an intermediate $\mathbf{M}$ point in the Brillouin zone. Here the
dispersion is linear along one direction and quadratic along the other. In the
presence of a Haldane flux, the system acquires topological properties, and
finite Chern numbers can be associated with the pairs of the conduction and
the valence bands. The valence band closer to the Fermi level ($E_{F}$)
possesses Chern numbers equal to $\pm 2$ and $\pm 1$, while the one further
away from $E_{F}$ corresponds to Chern numbers $\pm 1$. The conduction bands
are associated with similar properties, except their signs are reversed. The
Chern lobes shrink in the band-engineered model, and we find evidence of
multiple topological phase transitions, where the Chern numbers
discontinuously jump from $\pm 2$ to $\mp 2$, $\pm 1$ to $\mp 1$, $\pm 1$ to
$0$ to $\pm 2$ and $\pm 2$ to $\pm 1$. These transitions are supported by the
presence or absence of the chiral edge modes in a nanoribbon bilayer geometry
and the vanishing of the plateau in the anomalous Hall conductivity. Different
phases are further computed for different hopping amplitudes across the
layers, which shows shrinking of the Chern lobes for large interlayer
tunneling.

##  I Introduction

The Haldane model is a toy model, where it was shown that one can achieve
quantum Hall effect even in the absence of an external magnetic field in a
two-dimensional honeycomb lattice [1]. To achieve such a scenario, the time-
reversal symmetry (TRS) of the system needs to be broken which can be done via
chiral complex next nearest neighbour hopping amplitudes. The spectral bands
of such a system possess a non-zero topological invariant known as the Chern
number, and hence the system is known as a Chern insulator. Futhermore, the
band structure of a semi-infinite ribbon geometry hosts chiral edge modes
which is the signature of its topological character. Also, the system exhibits
quantum anomalous Hall effect which shows a plateau structure in the vicinity
of zero Fermi energy.

Haldane’s work has trigged an extensive study in both theoretical and
experimental fronts. for example, there have been reports of Haldane-like
spectrum and non-trivial phases realized in dice lattice [2], Kagomé lattice
[3, 4, 5, 6], checkerboard lattice [7], Lieb lattice [8, 9, 10, 11], buckled
lattice [12] etc. Experimentally, Haldane model has been realized in cold
atoms situated at the optical lattices where the complex second neighbour
hopping can be created by means of standing-wave laser beams [13, 14, 15],
ultracold fermions in the optical honeycomb lattices [16] etc. Also, a two-
dimensional honeycomb structure of Fe-based insulators, such as,
$X\mathrm{Fe_{2}(PO_{4})_{2}}$, with $X$ being either one of these, K, Cs, La
[17] demonstrate similar non-trivial topological phases with a non-zero Chern
number. Further, non-zero Chern numbers have also been found in acoustic Chern
insulators [18], the interface between the two trivial ferromagnetic
insulators EuO and GdN [19] etc.

In recent years, there have been studies of the Haldane model in coupled two-
dimensional systems, for example, the bilayer materials [20, 21, 22, 23]. In
parallel, there are studies on the band engineering in various systems, such
as, single layer graphene [24], spin Hall insulators [25] and a dice lattice
[26]. Such a band engineering has been incorporated via the introduction of an
anisotropy among the nearest neighbour (NN) hoppings. Such hopping
anisotropies have been included between the neighbouring sites lying along a
particular direction (say, $t_{1}$), while keeping the rest NN hoppings as $t$
in a honeycomb lattice. If the value of $t_{1}$ is varied, the band extrema
from the two Dirac points move closer to each other and they finally merge
with a vanishing band gap at the $\mathbf{M}$ point in the Brillouin zone (BZ)
for a particular value of $t_{1}$, namely, $t_{1}=2t$, which is called as the
semi-Dirac limit. During the process, the topological properties of the system
also vanishes at the gap closing hopping amplitude $t_{1}=2t$. It should be
noted that the band structure of the system in absence of the complex NNN
hopping (Haldane flux) shows semi-Dirac dispersion, that is, linear along the
$k_{x}$-direction and quadratic along the $k_{y}$-direction. Experimentally,
the semi-Dirac dispersion has been observed in many materials, such as,
multilayered structures of $\mathrm{TiO}_{2}/\mathrm{VO}_{2}$ [27, 28],
monolayer phosphorene in presence of doping and pressure [29, 30],
$\mathrm{BEDT}$-$\mathrm{TTF_{2}I_{3}}$ organic salts under pressure [31, 32],
black phosphorus doped with potassium atoms by means of in situ deposition
[33] etc. One can also achieve the semi-Dirac dispersion by applying an
uniaxial strain to a system which will change the bond length lying along a
particular direction which is parallel to the applied strain direction.
Therefore, the hopping energies along those directions get modified, while the
hopping along the other directions remains unaltered. Such method has been
employed in a monolayer honeycomb structure, such as, $\mathrm{Si_{2}O}$ which
yields a semi-Dirac dispersion [34].

However, the effect of band engineering in a multi-layered system, such as, a
bilayer graphene has never been studied. Needless to mention that bilayers
possess a richer phase diagram comprising of a larger parameter space. The
topological properties of such an engineered system are interesting since the
existence of the edge modes and the quantized Hall conductivity have never
been studied. A more interesting issue is that owing to larger number of bands
being present in the band structure of a bilayer system, higher values of
Chern numbers are realized. A higher Chern number implies higher value of the
anomalous Hall conductivity together with larger number of chiral edge modes
present in a semi-infinite system.

Higher Chern numbers are in general interesting and can be realized in a host
of systems, such as, in the Dirac [35] and semi-Dirac [36] systems in presence
of longer range hopping, multi-orbital triangular lattices [37], star lattices
or decorated honeycomb lattices [38], honeycomb lattices in presence of spin-
orbit coupling [39, 40], ultracold gases in triangular lattices [41, 42] etc.
Further, topological insulators doped with magnetic materials [43], Cr-doped
thin laminar sheets of $\mathrm{Bi}_{2}(\mathrm{Se,Te})_{3}$ [44] also
demonstrate higher values of the Chern numbers. Further,
$\mathrm{MnBi}_{2}\mathrm{Te}_{4}$ at high temperature [47, 46], multilayered
structure of doped (with magnetic materials) and undoped topological
insulators arranged alternatively [48], and in classical systems, such as,
sonic crystals prepared using acoustic components [45] show non-trivial phases
with higher Chern numbers.

(a)

(b)

Figure 1: A bilayer graphene is shown in (a) with the interlayer coupling
$t_{\perp}$ between the B sublattice of upper layer and the A sublattice of
lower layer. In both layers, the A and B sublattices are denoted by the red
and blue filled circles. In (b), the other planar hoppings are shown. To
properly see each sublattices in each layer, we have denoted the A and B
sublattices in lower layer with the circles in red and blue respectively. The
subscripts $l$ and $u$ in $\mathrm{A}_{l,u}$ and $\mathrm{B}_{l,u}$ refers to
lower and upper layer respectively. All the bondings and NNN hoppings in the
lower layers are shown by the dashed lines and dashed arrows respectively. The
NN hopping strength along the $\boldsymbol{\delta}_{1}$ direction (shown via
the yellow arrow) is $t_{1}$, while it is $t$ along the
$\boldsymbol{\delta}_{2,3}$ directions ($\boldsymbol{\delta}_{i}$ are defined
in text). The NNN hopping is $t_{2}e^{i\phi}$ ($t_{2}e^{-i\phi}$) for the
clockwise (anti-clockwise) direction.

In this work, we focus on a bilayer graphene with broken TRS, that is, a
coupled bilayer Haldane model. The stacking of the two layers is assumed in
such a way that the B subllatice of the upper layer lies exactly above the A
sublattice of the lower layer. Such stacking is known as the AB stacking or
the Bernal stacking. We shall see that the Chern numbers associated with
various bands reveal interesting properties. For example, some of the bands
possess both Chern numbers $\pm 2$ and $\pm 1$, while the rests are associated
with Chern numbers $\pm 1$. Such a scenario needs to be assessed for a band
engineered system. Specifically, we wish to address the ramifications of the
band deformation caused via asymmetric hopping amplitudes on the topological
properties and ascertain whether such deformation induces a topological phase
transition. In our bilayer model, the band engineering is incorporated via
asymmetric NN hopping amplitudes in each of the layers, while the tunneling
amplitude across the layers is left unaltered.

Our subsequent discussions have been arranged as follows. Sec. II introduces
the tight binding Hamiltonian of a bilayer graphene. Sec. III discusses the
band structure of the system with the interlayer coupling ($t_{\perp}$) and
the anisotropic NN hopping amplitudes ($t_{1}$) as parameters. Sec. IV deals
with the phase diagrams that are obtained by computing the Chern numbers
associated with the bands. In Sec. V, the presence (or absence) of the chiral
edge modes in a ribbon geometry are presented. Next, the numerical
computations of the anomalous Hall conductivity are shown in Sec. VI. Finally,
a brief summary of the results are included in the concluding section (Sec.
VII).

##  II The Hamiltonian

(a) (b) (c) (d) (e) (f) (g) (h) (i)

Figure 2: The band structure in absence of $t_{2}$ $(t_{2}=0)$ is shown along
the $k_{x}$-axis (at $k_{y}a_{0}=2\pi/3$) for (a) $t_{1}=t$, (b) $t_{1}=1.5t$,
(c) $t_{1}=2t$, and (d) $t_{1}=2.2t$. Similarly, the dispersions in presence
of $t_{2}$ $(t_{2}=0.1t)$ are depicted for (e) $t_{1}=t$, (f) $t_{1}=1.5t$,
(g) $t_{1}=1.8t$, (h) $t_{1}=2t$, and (i) $t_{1}=2.2t$. The values of the
other parameters are $t_{\perp}=0.5t$ and $\phi_{l}=\phi_{u}=\pi/2$.

(a) (b) (c) (d) (e) (f) (g) (h) (i)

Figure 3: The band structure for $t_{2}=0$ is shown along the $k_{x}$-axis (at
$k_{y}a_{0}=2\pi/3$) for (a) $t_{1}=t$, (b) $t_{1}=1.5t$, (c) $t_{1}=2t$, and
(d) $t_{1}=2.2t$. While, the spectra for a non-zero $t_{2}$ $(t_{2}=0.1t)$ are
depicted for (e) $t_{1}=t$, (f) $t_{1}=1.5t$, (g) $t_{1}=1.8t$, (h)
$t_{1}=2t$, and (i) $t_{1}=2.2t$. The values of $t_{\perp}$, $\phi_{l}$ and
$\phi_{u}$ are fixed at $0.1t$, $\pi/2$ and $\pi/2$ respectively.

A tight-binding Hamiltonian of a bilayer honeycomb lattice can be written as
follows,

$\displaystyle H=$ $\displaystyle\sum_{p\in l,u}\left[\sum_{\langle
ij\rangle}t_{ij}c_{i}^{p\dagger}c^{p}_{j}+t_{2}\sum_{\langle\langle
im\rangle\rangle}e^{i\phi^{im}_{p}}c_{i}^{p\dagger}c^{p}_{m}+\mathrm{h.c.}\right]$
$\displaystyle+\left[t_{\perp}\sum_{{\left\langle
q,r\right\rangle}_{\perp}}c_{q}^{l\dagger}c^{u}_{r}+\mathrm{h.c.}\right]$ (1)

where $c_{i}^{p\dagger}(c_{i}^{p})$ is the creation (annihilation) operator
corresponding to site $i$ which belongs to the layer $p$. Here $p=l,u$
represent the lower and the upper layers respectively. The first term in the
right hand side denotes the nearest neighbour (NN) hopping with the amplitude
$t_{ij}$ being either $t_{1}$ when $i$ and $j$ sites lie along the
$\boldsymbol{\delta}_{1}=a_{0}(0,1)$ direction, or $t$ when they lie along the
$\boldsymbol{\delta}_{2}=a_{0}(\sqrt{3}/2,-1/2)$ and
$\boldsymbol{\delta}_{3}=a_{0}(\sqrt{3}/2,-1/2)$ directions as shown in Fig.
1. The second term represents the complex next nearest neighbour (NNN) hopping
with the amplitude $t_{2}$ and a phase $\phi^{im}_{l,p}$. We have labeled the
Haldane flux corresponding to the lower and upper layers as $\phi^{im}_{l}$
and $\phi^{im}_{u}$ respectively. If an electron hops in the counter-clockwise
direction, $\phi^{im}_{l,p}$ assumes a positive sign, while for the clockwise
direction, it acquires negative sign. The third term is the hopping between
the two layers with the coupling strength $t_{\perp}$. It should be kept in
mind that the interlayer hopping is between the B sublattice on layer $u$
($r\in\mathrm{B}_{u}$) and the A sublattice on layer $l$
($q\in\mathrm{A}_{l}$) (AB or Bernal stacking). In our calculations, we have
varied $t_{1}$ in both the layers from a value $t$ to $2t$ (semi-Dirac) and
even considered $t_{1}>2t$.

Now, we Fourier transform the Hamiltonian and write them in the four
sublattice basis, namely,
$\\{\mathrm{A}_{l},\mathrm{B}_{l},\mathrm{A}_{u},\mathrm{B}_{u}\\}$ in the
following way,

$H(\mathbf{k})=\begin{pmatrix}h^{+}_{z}(\mathbf{k},\phi_{l})&h_{xy}(\mathbf{k},t_{1})&0&t_{\perp}\\\
h_{xy}^{*}(\mathbf{k},t_{1})&h_{z}^{-}(\mathbf{k},\phi_{l})&0&0\\\
0&0&h^{+}_{z}(\mathbf{k},\phi_{u})&h_{xy}(\mathbf{k},t_{1})\\\
t_{\perp}&0&h_{xy}^{*}(\mathbf{k},t_{1})&h^{-}_{z}(\mathbf{k},\phi_{u})\\\
\end{pmatrix}$ (2)

where $h_{z}^{\pm}$ are defined as,
$h^{+}_{z}(\mathbf{k},\phi_{p})=h_{0}(\mathbf{k},\phi_{p})\pm
h_{z}(\mathbf{k},\phi_{p}$). The element $h_{xy}(\mathbf{k},t_{1})$ has the
following form,
$h_{xy}(\mathbf{k},t_{1})=h_{x}(\mathbf{k},t_{1})-ih_{y}(\mathbf{k},t_{1})$.
The expressions for the $h_{i}$s can be written as,

$\displaystyle
h_{0}(\mathbf{k},\phi_{p})=2t_{2}\cos\phi_{p}\left\\{2\cos\frac{\sqrt{3}k_{x}}{2}\cos\frac{3k_{y}}{2}+\cos\sqrt{3}k_{x}\right\\}$
(3) $\displaystyle
h_{z}(\mathbf{k},\phi_{p})=-2t_{2}\sin\phi_{p}\left\\{2\sin\frac{\sqrt{3}k_{x}}{2}\cos\frac{3k_{y}}{2}-\sin\sqrt{3}k_{x}\right\\}$
(4) $h_{x}(\mathbf{k},t_{1})=\left\\{t_{1}\cos
k_{y}+2t\cos\frac{k_{y}}{2}\cos\frac{\sqrt{3}k_{x}}{2}\right\\},$ (5)

and

$h_{y}(\mathbf{k},\tilde{t})=\left\\{-t_{1}\sin
k_{y}+2t\sin\frac{k_{y}}{2}\cos\frac{\sqrt{3}k_{x}}{2}\right\\},$ (6)

Throughout our work, the amplitude of the NNN hopping $t_{2}$ is kept fixed at
$0.1t$, and two different values of the interlayer hopping strength are
chosen, namely, $t_{\perp}=0.5t$ and $t_{\perp}=0.1t$ [49]. The values of
$\phi_{l}$ and $\phi_{u}$ are taken such that $\phi_{l}=\phi_{u}=\pi/2$. Now,
for $\phi_{u}=\phi_{l}$, we obtain the following dispersion relation,

$\displaystyle
E^{c}_{\pm}=\left[h_{0}+\sqrt{\frac{t_{\perp}^{2}}{2}+|h_{xy}|^{2}+h_{z}^{2}\pm\frac{t_{\perp}}{2}\sqrt{t_{\perp}^{2}+4h_{xy}^{2}}}\right]$
(7) $\displaystyle
E^{v}_{\pm}=\left[h_{0}-\sqrt{\frac{t_{\perp}^{2}}{2}+|h_{xy}|^{2}+h_{z}^{2}\pm\frac{t_{\perp}}{2}\sqrt{t_{\perp}^{2}+4h_{xy}^{2}}}\right]$
(8)

where $E^{c}_{\pm}$ denote the two conduction bands and $E^{v}_{\pm}$ are the
two valence bands for a bilayer.

(a) (b) (c) (d)

Figure 4: The phase diagrams corresponding to the lowest occupied band, that
is, band-v1 is presented for $t_{\perp}=0.5t$. The white regions denote the
trivial phase with Chern number as zero, while the colored regions indicate
the non-trivial phase with the non-zero Chern numbers. The non-zero values are
indicated at the top of the figure.

(a) (b) (c) (d)

Figure 5: The phase diagrams corresponding to band-v2 is presented for
$t_{\perp}=0.5t$. The white regions denote the trivial phases with zero Chern
number, while the colored regions indicate the non-trivial phase with the non-
zero Chern numbers. Again the values are indicated at the top of the figure.

##  III Spectral properties

In this section, we discuss how the spectral properties evolve as we
interpolate between the Dirac and the semi-Dirac limits. We show the band
structure for two different values of $t_{\perp}$. The first one is for
$t_{\perp}=0.5t$ as shown in Fig. 2. As can be seen, there are four bands
which we have labeled as follows. The upper conduction band is labeled as
band-c1, while the lower conduction band is band-c2. Similarly, the lower and
the upper valence bands are labeled as band-v1 and band-v2 respectively. When
$t_{2}=0$ (no Haldane flux), band-c2 and band-v2 touch each other at the Fermi
level at the $\mathbf{K}$ and $\mathbf{K^{\prime}}$ points (see Figs.
2(a)-2(d)). These points are referred to as the Dirac points. Further, with
the increase in the value of $t_{1}$, we deviate from the Dirac limit, and the
band touching points move close to each other which finally merge at
$t_{1}=2t$. Beyond this value, that is, for $t_{1}>2t$, a gap opens up at the
$\mathbf{M}$ point. Now, if we switch on $t_{2}$ (see Figs. 2(e)-2(i)), the
spectral gap remains open for $t\leq t_{1}<2t$ and $t_{1}>2t$, while the gap
vanishes exactly at the semi-Dirac limit, namely, $t_{1}=2t$. The gap closing
scenario of the bilayer graphene is thus similar to the case of single layer
graphene, where the energy gap between the conduction and valence band
vanishes at the semi-Dirac limit, that is, at $t_{1}=2t$ [24].

Further, we have presented band structure in Fig. 3 for a smaller value of
$t_{\perp}$, namely, $t_{\perp}=0.1t$. It is obvious from Eqs. 7 and 8 that
the separation among the conduction bands (band-c1 and band-c2) and that among
the valence bands (band-v1 and band-v2) decreases with decrease in
$t_{\perp}$. Moreover, the low energy dispersions of band-c2 and band-v2 about
the band touching points have a linear behaviour which was quadratic for
$t_{\perp}=0.5t$. Thus, the massive electrons become progressively massless as
we lower the value of $t_{\perp}$. Further, with the decrease in $t_{\perp}$,
the spectral gap between band-c2 and band-v2 increases. For example, when
$t_{\perp}=0.1t$ the band gap is $\Delta E_{g}\simeq 1.0390t$ and $0.3124t$
for $t_{1}=t$ and $1.8t$ respectively. While for $t_{\perp}=0.5t$, $\Delta
E_{g}\simeq 1.0335t$ and $0.1406t$ for $t_{1}=t$ and $1.8t$ respectively. Thus
the difference in energy is more noticeable as we move towards the semi-Dirac
limit, that is, at large values of $t_{1}$.

##  IV Chern number and phase diagram

(a) (b) (c) (d)

Figure 6: The phase diagrams corresponding to band-v2 are shown in (a) and
(b), and those for band-v1 are presented in (c) and (d). In (a) and (c), the
$\eta$ points are used and shown along the $L_{1}$ and $L_{2}$ lines, whereas
the $\gamma$ and $\chi$ points are marked along the $L_{3}$ and $L_{4}$ lines
in both (b) and (d). Along those lines multiple phase transitions occur. For
example, along $L_{3}$, the Chern number corresponding to band-v2 has values
$+1$, $0$, $+2$, $0$, $+1$ at the points $\gamma_{1}$, $\gamma_{3}$,
$\gamma_{5}$, $\gamma_{7}$ and $\gamma_{9}$ respectively. The phase
transitions take place at $\gamma_{2}$, $\gamma_{4}$, $\gamma_{6}$ and
$\gamma_{8}$, where band-v2 touches either band-v1 or band-c2. The values of
$t_{\perp}$ and $t_{1}$ are taken as $0.5t$ and $t_{1}$ respectively.

In this section, we calculate the Chern number as a function of the Haldane
flux of the two layers. Owing to the broken TRS, the bands possess non-zero
Chern numbers, which can be calculated by integrating the Berry curvature over
the BZ [50, 51].

$\displaystyle C$ $\displaystyle=$
$\displaystyle\frac{1}{2\pi}\int\int_{\mathrm{BZ}}\Omega(k_{x},k_{y})\mathrm{d}k_{x}\mathrm{d}k_{y}$
(9)

where $\Omega(k_{x},k_{y})$ is the $z$-component of the Berry curvature [52],
which is obtained from the following relation.

$\displaystyle\Omega(k_{x},k_{y})=-2i\mathrm{Im}\left[\left<\frac{\partial\psi(k_{x},k_{y})}{\partial
k_{x}}\right.\left|\frac{\partial\psi(k_{x},k_{y})}{\partial
k_{y}}\right>\right]$ (10)

where $\psi(k_{x},k_{y})$ is the periodic part of the Bloch wave corresponding
to the Hamiltonian defined in Eq. 2, and $\mathrm{Im}$ denotes the imaginary
part. Hence, we calculate the Chern numbers as a function of the fluxes
$\phi_{l}$ and $\phi_{u}$ corresponding to the lower and the upper layers
respectively for various values of $t_{1}$ as shown in Fig. 4. Here, the value
of $t_{\perp}$ is chosen to be $0.5t$ and the phase diagrams shown correspond
to band-v1. We have denoted the Chern insulating regions by two colors. The
regions in red denote $C=+1$ phase, while the blue ones denote $C=-1$ phases.
The trivial phases with $C=0$ are shown by the white regions. It is evident
from Fig. 4(a) that the areas of the Chern insulating regions are maximum for
$t_{1}=t$ (Dirac case). An engineering of the band structure, that is, with
the increase in the value of $t_{1}$, the area of the topological regions
(called as the Chern lobes) gradually shrink. We have shown the phase diagram
till a certain value, namely, $t_{1}=1.9t$ (see Fig. 4(d)), beyond which the
topological regions can hardly be seen. When $t_{1}$ becomes equal to $2t$,
the Chern number ($C$) vanishes completely for all values of $\phi_{l}$ and
$\phi_{u}$ owing to a gapless scenario between band-c2 and band-v2. Although
band-v1 remains separated from the band-v2, the Chern number still vanishes.
For $t_{1}>2t$, a gap opens up, however, the Chern numbers continue to be
zero, and thus the gap is trivial.

Further, we have presented the phase diagrams corresponding to band-v2 (the
one closer to the Fermi level) in Fig. 5. As can be seen, additional phases
with higher Chern number ($C=\pm 2$) appear. We have denoted the $C=+2$ and
$C=-2$ phase with cyan and green colors respectively. The red and blue colors
continue to denote $C=+1$ and $C=-1$ phases respectively. Thus, both $C=\pm 2$
and $C=\pm 1$ phases occur at different parameter values in the same phase
diagram. Further, the topological regions shrinks with the increase in
$t_{1}$, and finally vanishes at $t_{1}=2t$, where the gap between the band-v2
and band-c2 vanishes. For $t_{1}>2t$, the gap reopens, but the Chern number
remains zero for all values of $\phi_{l}$ and $\phi_{u}$. The phase diagrams
for band-c1 and band-c2 are identical in shape to those of band-v1 and band-v2
respectively, except the Chern numbers have opposite signs.

In order to visualize the gap closing scenario corresponding to different
phase transitions occurring in the phase diagrams, the band structures are
presented in Fig. 7 for a particular value of $t_{1}$ and $t_{\perp}$, namely,
$t_{1}=t$ and $t_{\perp}=0.5t$. The values of $\phi_{l}$ and $\phi_{u}$ are
such that they lie along the four lines, namely, $L_{1}$, $L_{2}$, $L_{3}$ and
$L_{4}$ in the phase diagrams depicted in Fig. 6, and are denoted by
$\eta_{i}\,(i=1,\dots,6)$, $\gamma_{j}\,(j=1,\dots,9)$, and
$\chi_{s}\,(s=1,2,3)$. Along $L_{1}$, a topological phase transition occur
between $C=+2$ and $C=-2$ corresponding to band-v2, while the transition
between $C=+1$ and $C=-1$ occur along $L_{2}$ for both band-v2 and band-v1.
These results have to be understood in conjunction with the corresponding band
structures as shown in Figs. 7(a)-7(c) and 7(d)-7(f) respectively. The band
structures corresponding to $\eta_{1}$ and $\eta_{3}$ points are identical and
the Chern numbers corresponding to band-v2 are $+2$ and $-2$ respectively
(Fig. 6(a)). At $\eta_{2}$, band-v2 and band-c2 touch each other at both the
Dirac points (Fig. 7(b)), and hence there is a phase transition at $\eta_{2}$.
However, band-v1 remains isolated from band-v2 at these $\eta$ points, and the
Chern numbers are zero along $L_{1}$ as evident from its phase diagram (Fig.
6(c)). Further, the band structures corresponding to $\eta_{4}$ and $\eta_{6}$
have similar features, however in this case, $C$ has values $+1$ and $-1$
respectively corresponding to band-v2, while for band-v1, $C$ has the same
magnitude, but are of opposite signs. At the phase transition occurring at
$\eta_{5}$, band-v2 and band-v1 touch each other at the $\mathbf{K}$ point in
the BZ (Fig. 7(e)), and hence for both the bands, a topological phase
transition takes place at this point.

Further, along $L_{3}$, again multiple phase transitions occur (see Figs. 6(b)
and 6(d)), and the corresponding dispersions are shown in Fig. 7(g)-7(o). At
$\gamma_{1}$, band-v2 and band-v1 show $C=+1$ and $C=-1$ respectively which
drops to zero at $\gamma_{2}$ and hence the gap between those bands close at
the $\mathbf{K}$ point as shown in Fig. 7(h). At $\gamma_{3}$, the gap
reopens, but the Chern numbers corresponding to these bands remain zero. The
gap between band-v2 and band-c2 vanishes at $\gamma_{4}$ where again a phase
transition takes place, since along the line connecting $\gamma_{4}$ and
$\gamma_{6}$, the Chern number has a value $+2$. The band structure at an
intermediate point, namely, $\gamma_{5}$ has been shown in Fig. 7(k).
Similarly, phase transitions take place at $\gamma_{6}$ and $\gamma_{8}$,
where the gaps vanish at the $\mathbf{K^{\prime}}$ point. At $\gamma_{7}$ and
$\gamma_{9}$, $C$ assumes values zero and $+1$ respectively corresponding to
band-v2. It should be noted that band-v1 shows vanishing of the Chern number
between $\gamma_{2}$ and $\gamma_{8}$ segments (see fig. 6(d)) and hence it
never touches band-v2 which results in absence of any phase transition.

Now, we show the phase transitions between $C=-2$ and $C=-1$ phase along
$L_{4}$. The corresponding band structures are shown in Fig. 7(p)-7(r). At
$\chi_{2}$, band-v2 and band-v1 remain isolated from each other, however, they
possess Chern numbers $C=-2$ and $C=0$ respectively. At $\chi_{1}$ and
$\chi_{3}$, these two bands touch each other at the $\mathbf{K^{\prime}}$ and
the $\mathbf{K}$ points in the BZ respectively, where topological phase
transitions take place. Beyond $\chi_{1}$ and $\chi_{3}$, the gap reopens and
both the bands possess non-trivial phases with $C=-1$.

Further, along the $\phi_{u}=-\phi_{l}$ line, a semi-metallic phase exists for
all the bands. In the vicinity of $\phi_{u}=\phi_{l}$, only the phase diagrams
of band-v1 show trivial regions with $C=0$, however, those for band-v2
demonstrate non-trivial phases either with $C=+2$ or $C=-2$.

Moreover, in order to see the effects of $t_{\perp}$ on the topological
phases, we have shown the phase diagrams corresponding to band-v1 and band-v2
in Figs. 8(a)-8(d) and 8(e)-8(h) respectively. It is evident that the areas of
Chern insulating regions are enhanced corresponding to lower values of
$t_{\perp}$. Also, the shape of the topological regions are different from
those for the $t_{\perp}=0.5t$ case. Further, the areas of $C=\pm 2$ regions
in the phase diagram corresponding to band-v2 are mostly spanned by $C=\pm 1$
regions. However, the feature that remains unaltered is the trivial phase
along $\phi_{u}=\pm\phi_{l}$ lines for band-v1 and $\phi_{u}=\phi_{l}$ line
for band-v2. For both $t_{\perp}=0.5t$ and $t_{\perp}=0.1t$, the Chern
insulating regions gradually shrink with the increase in the value of $t_{1}$
and finally vanish at the semi-Dirac limit, namely, $t_{1}=2t$.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

(j)

(k)

(l)

(m)

(n)

(o)

(p)

(q)

(r)

Figure 7: The band structures corresponding to the points
$\eta_{1}$-$\eta_{6}$ (shown in Figs. 6(a) and 6(c)) are depicted in (c)-(h).
The spectra for the points $\gamma_{1}$-$\gamma_{9}$ (shown in Figs. 6(a) and
6(c)) are shown in (i)-(q), and for the points $\chi_{1}$-$\chi_{3}$ are
presented in (p)-(r). The values of $t_{1}$ and $t_{\perp}$ for all the band
structures are kept fixed at $t$ and $0.5t$ respectively.

(a) (b) (c) (d) (e) (f) (g) (h)

Figure 8: The phase diagrams corresponding to band-v1 are shown in (a)-(d),
while those for band-v2 are shown in (e)-(h). The value of $t_{\perp}$ is kept
fixed at $0.1t$. The values of $t_{1}$ are such that $t_{1}=t$ in (a) and (e),
$t_{1}=1.5t$ in (b) and (f), $t_{1}=1.8t$ in (c) and (g), and $t_{1}=1.9t$ in
(d) and (h). The white regions denote trivial phases with zero Chern numbers,
while the colored regions indicate the non-trivial phases with the non-zero
Chern numbers. The values are indicated at the top of the figure.

##  V Edge states

(a)

(b)

(c)

(d)

Figure 9: The edge state spectra are shown for (a) $t_{1}=t$, (b)
$t_{1}=1.5t$, (c) $t_{1}=2t$, and (d) $t_{1}=2.2t$. The green shaded regions
represent the bulk gap in (a), (b) and (d) (np bulk gap in (c)). The Fermi
levels ($E_{F}$) are denoted by the red dashed lines, which are shown to be
present in the bulk gap. $E_{F}$ intersects the edge modes at the points
denoted by the green dots as shown in (a) and (b). For these, the edge
currents are shown by the green arrows in the yellow panels located at the top
right corner, which represent parts of the semi-infinite ribbon.

To show the existence (and their vanishing) of the edge modes, in this
section, we show the band structure of the system for semi-infinite
nanoribbon. The ribbon has a finite width along the $y$-direction, while it is
infinite along the $x$-direction [54, 53]. Further, we label the sites along
the $y$-direction as A${}_{1}^{l}$, B${}_{1}^{l}$, A${}_{2}^{l}$,
B${}_{2}^{l}$, …. A${}_{N}^{l}$, B${}_{N}^{l}$, A${}_{1}^{u}$, B${}_{1}^{u}$,
A${}_{2}^{u}$, B${}_{2}^{u}$, …. A${}_{N}^{u}$, B${}_{N}^{u}$. Since, the
periodicity along the $x$-direction remains preserved, we can Fourier
transform the operators along that direction. This results in set of four
coupled equations as shown below.

$\displaystyle E_{k}a_{k,n}^{u}=$
$\displaystyle-\left[t\left\\{1+e^{(-1)^{n}ik}\right\\}b_{k,n}^{u}+t_{1}b_{k,n-1}^{u}\right]$
(11)
$\displaystyle-2t_{2}\left[\cos(k+\phi)a_{k,n}^{u}+e^{(-1)^{n}\frac{ik}{2}}\times\right.$
$\displaystyle\left.\cos\left(\frac{k}{2}-\phi\right)\\{a_{k,n-1}^{u}+a_{k,n+1}^{u}\\}\right]$

$\displaystyle
E_{k}b_{k,n}^{u}=-\left[t\left\\{1+e^{(-1)^{n+1}ik}\right\\}a_{k,n}^{u}+t_{1}a_{k,n+1}^{u}\right]$
(12)
$\displaystyle-2t_{2}\left[e^{(-1)^{n+1}\frac{ik}{2}}\times\cos\left(\frac{k}{2}+\phi\right)\\{a_{k,n-1}^{u}+a_{k,n+1}^{u}\\}+\right.$
$\displaystyle\left.\cos(k-\phi)b_{k,n}^{u}\right]+t_{\perp}\left[\xi_{1}e^{-ik}+\xi_{2}\right]a_{n}^{l}$
$\displaystyle
E_{k}a_{k,n}^{l}=-\left[t\left\\{1+e^{(-1)^{n}ik}\right\\}b_{k,n}^{l}+t_{1}b_{k,n-1}^{l}\right]$
(13)
$\displaystyle-2t_{2}\left[e^{(-1)^{n}\frac{ik}{2}}\times\cos\left(\frac{k}{2}-\phi\right)\\{a_{k,n-1}^{l}+a_{k,n+1}^{l}\\}\right.$
$\displaystyle\left.\cos(k+\phi)a_{k,n}^{l}+\right]+t_{\perp}\left[\xi_{1}e^{ik}+\xi_{2}\right]b_{n}^{u}$
$\displaystyle E_{k}b_{k,n}^{l}=$
$\displaystyle-\left[t\left\\{1+e^{(-1)^{n+1}ik}\right\\}a_{k,n}^{l}+t_{1}a_{k,n+1}^{l}\right]$
(14)
$\displaystyle-2t_{2}\left[\cos(k-\phi)b_{k,n}^{l}+e^{(-1)^{n+1}\frac{ik}{2}}\times\right.$
$\displaystyle\left.\cos\left(\frac{k}{2}+\phi\right)\\{a_{k,n-1}^{l}+a_{k,n+1}^{l}\\}\right]$

where $a_{k,n}^{l,u}$ and $b_{k,n}^{l,u}$ are the amplitudes of the wave
functions corresponding to the sublattices A and B respectively. The
superscripts $l$ and $u$ refers to lower and upper layers respectively. Here
$k=\sqrt{3}k_{x}a_{0}$ is the dimensionless momentum and $n$ denotes the site
index which assumes integer values in the range $[1:N]$ with $N$ being the
total number of unit cells along the $y$-direction. We chose $N$ as 128, which
gives the width to be $79\sqrt{3}a_{0}$. In Eqs. 12 and 13, $\xi_{1}$ and
$\xi_{2}$ denote quantities that depend on the site index, $n$ via
$\xi_{1}=[1-(-1)^{n}]/2$ and $\xi_{2}=[1+(-1)^{n}]/2$ respectively.

By solving Eqs. 11, 12, 13, and 14, we obtain the band structure of the
nanoribbon for various values of $t_{1}$ as presented in Fig. 9. It can be
noticed that a pair of edge modes from the valence bands (band-v2) traverse
the Fermi level, $E_{F}$ (shown via the red dashed line) and merge with the
conduction bands (band-c2) and another pair crosses the Fermi level in the
opposite direction. Such crossing of the edge modes leads to a quantized Hall
conductivity should the Fermi level lies in the bulk gap. $E_{F}$ intersects
the edge modes (see Figs. 9(a) and 9(b)) at four points (marked by the green
dots), whose corresponding edge currents are shown by the green arrows in the
yellow panels located at the top right corner of the plots. The yellow panels
represent a part of the semi-infinite ribbon. Since the velocity of electrons
are proportional to the slope of the band structure, that is, $\partial
E/\partial k$, there exists a pair of edge currents at each edge that moves in
the same direction. However, such pairs of currents propagate in opposite
directions at the two edges of the ribbon. Hence, these modes are called
chiral edge modes.

It should be noted that because of a pair of chiral edge modes, we should
obtain the Hall conductivity quantized with a plateau at $2e^{2}/h$, with the
factor ‘2’ arise due to doubling of the number of chiral edge modes [55]. Such
chiral edge modes exist as long as the value of $t_{1}$ remains lesser than
$2t$. Since the bulk gap vanishes at $t_{1}=2t$ (see Fig. 9(c)), the edge
current vanishes. For $t_{1}>2t$, the edge modes get detached from the bulk
bands as shown in Fig. 9(d) for $t_{1}=2.2t$, thereby resulting in a zero edge
current. These results are consistent with the corresponding Chern numbers
obtained in the phase diagram. For example, we observe the non-zero edge
currents for $t_{1}<2t$ and the corresponding Chern number is found to be
$C=|2|$. For $t_{1}>2t$, the Chern numbers vanish, and so the edge currents.
The figures presented here are for $t_{\perp}=0.5t$. For $t_{\perp}=0.1t$, we
observe similar features in the spectrum, except that the bulk gaps get
reduced. We have skipped the discussion of the latter for brevity.

##  VI Hall conductivity

(a)

(b)

Figure 10: The anomalous Hall conductivities are shown as a function of
$E_{F}$ for various values of $t_{1}$ in (a) and (b) for $t_{\perp}=0.5t$ and
$t_{\perp}=0.1t$ respectively. The plateau width decreases as $t_{1}$ deviates
from $t$.

In this section, we calculate the anomalous Hall conductivity as function of
the Fermi energy $E_{F}$. The prerequisite is the computations of the Berry
curvature using Eq. 10 and hence use the following formula to calculate the
anomalous Hall conductivity ($\sigma_{xy}$) [56, 57], namely,

$\sigma_{xy}=\frac{\sigma_{0}}{2\pi}\sum_{\lambda}\int\frac{\mathrm{d}k_{x}\mathrm{d}k_{y}}{(2\pi)^{2}}f\left(E^{\lambda}_{k_{x},k_{y}}\right)\Omega(k_{x},k_{y})$
(15)

where
$f(E^{\lambda}_{k_{x},k_{y}})=\left[1+\exp{\\{(E^{\lambda}_{k_{x},k_{y}}-E_{F})/K_{B}T\\}}\right]^{-1}$
is the Fermi-Dirac distribution function at an energy
$E^{\lambda}_{k_{x},k_{y}}$. Here $E_{F}$ refers to the Fermi energy and $T$
is the absolute temperature. The energy is denoted by
$E^{\lambda}_{k_{x},k_{y}}$ with $\lambda$ being the band index. The constant
term $\sigma_{0}$ is equal to $e^{2}/h$ which sets the scale for
$\sigma_{xy}$. Now, we compute $\sigma_{xy}$ numerically as a function of
$E_{F}$ at zero temperature ($T=0$) for various values of $t_{1}$ as shown in
Fig. 10.

As can be seen from Fig. 10(a), when the Fermi energy, $E_{F}$ lies in the
bulk gap, $\sigma_{xy}$ becomes quantized at a value $2\sigma_{0}$. The width
of the plateau is equal to the width of the bulk gap in the dispersion
spectrum of Fig. 2. As soon as $E_{F}$ intersects the bands (either both the
conduction or both the valence bands), $\sigma_{xy}$ starts to decrease since
the integral is performed over the occupied states. This also results in
diminishing of the plateau width with increase in the value of $t_{1}$. This
happens because the energy gap between the band-c2 and the band-v2 shrinks.
The plateau and the Hall conductivity vanish completely at $t_{1}=2t$, where
the spectrum becomes gapless. For the hopping asymmetry engineered beyond the
semi-Dirac limit, that is, $t_{1}>2t$, the bands become gapped again, however,
the Hall conductivity remains zero. These results are consistent with their
corresponding information coming from the Chern numbers. The Hall plateaus are
observed as long as the system remains a Chern insulator, that is, for
$t_{1}<2t$. Further, the factor ‘2’ in ($2\sigma_{0}$) denotes the value of
Chern number (and also the edge modes) which vanishes for $t_{1}>2t$.

We have also presented the Hall conductivity for a smaller value of
$t_{\perp}$, namely, $t_{\perp}=0.1t$ in Fig. 10(b). In this case, the plateau
widths corresponding to different values of $t_{1}$ are larger as compared to
that for the $t_{\perp}=0.5t$ case, since the corresponding band gaps are
enhanced as shown in Fig. 3. However, similar to the previous case, the
plateau width decreases with the increase of $t_{1}$, which finally vanishes
at $t_{1}=2t$ and beyond. Thus, a topological phase transition takes place
across the gap closing point at the semi-Dirac limit, namely, $t_{1}=2t$.

##  VII Conclusion

We have investigated the topological properties for a band engineered bilayer
Haldane. By tuning one of the three NN hopping amplitudes, the band extrema,
which were located at the $\mathbf{K}$ and the $\mathbf{K^{\prime}}$ points
for the Dirac case, migrate towards each other and finally merge at an
intermediate $\mathbf{M}$ point in the BZ in the semi-Dirac limit, that is, at
$t_{1}=2t$. We have calculated the Chern numbers for various values of $t_{1}$
and plotted them in the $\phi_{u}$-$\phi_{l}$ plane which demonstrates that
the higher Chern numbers ($C=\pm 2$) are associated only with band-v2.
However, the Chern numbers corresponding to both both bands vanish, that is,
there are topological phase transitions, where the Chern numbers
discontinuously change from $C=\pm 2$ to $C=0$ and $C=\pm 1$ to $C=0$, across
the semi-Dirac point $t_{1}=2t$. Also, there are multiple phase transitions in
the phase diagram, such as, $+2\rightarrow-2$, $+1\rightarrow-1$, $\pm
2\rightarrow\pm 1$ and $\pm 2\rightarrow 0\rightarrow 1$. These phase
transitions are confirmed by the opening and closing of the energy gaps (semi-
metallic phase) in the dispersion spectrum. Further, we have also computed the
band structure of a nanoribbon, where we observe a pair of chiral edge modes
along the edges of the ribbon exist as long as $t_{1}$ remains lesser than
$2t$. Also for the anomalous Hall conductivity, the width of the quantized
plateau at $2\sigma_{0}$ gradually decreases with increase in $t_{1}$, which
finally vanishes at $t_{1}=2t$. Thus, a bilayer Haldane model, similar to its
monolayer analogue, exhibits a topological phase transition at the semi-Dirac
point. However, here we have larger values of the Chern number and doubling of
the edge modes at the edges of the bilayer nanoribbon. Further, the phase
transitions are supported by the vanishing of Chern number, chiral edge modes
and the anomalous Hall conductivity.

## References

  * [1] F. D. M. Haldane, Phys. Rev. Lett. 61, 2015 (1988).
  * [2] B. Dey, P. Kapri, O. Pal, and T. K. Ghosh, Phys. Rev. B 101, 235406 (2020).
  * [3] K Ohgushi, S. Murakami, and N. Nagaosa, Phys. Rev. B 62, R6065(R) (2000).
  * [4] Y. Xiao, V. Pelletier, P. M. Chaikin, and D. A. Huse, Phys. Rev. B 67, 104505 (2003).
  * [5] H. -M. Guo and M. Franz, Phys. Rev. B 80, 113102 (2009).
  * [6] X. -P. Liu, W. -C. Chen, Y. -F. Wang and C. -D. Gong J. Phys. Condens. Matter 25, 305602 (2013).
  * [7] K. Sun, H. Yao, E. Fradkin, and S. A. Kivelson Phys. Rev. Lett. 103, 046811 (2009).
  * [8] C. Weeks and M. Franz, Phys. Rev. B 82, 085310 (2010).
  * [9] V. Apaja, M. Hyrkäs, and M. Manninen, Phys. Rev. A 82, 041402(R) (2010).
  * [10] N. Goldman, D. F. Urban, and D. Bercioux, Phys. Rev. A 83, 063601 (2011).
  * [11] W. -F. Tsai, C. Fang, H. Yao, and J. Hu, New J. Phys. 17, 055016 (2015).
  * [12] A. Wright, Sci. Rep. 3, 2736 (2013).
  * [13] L. B. Shao, S. -L. Zhu, L. Sheng, D. Y. Xing, and Z. D. Wang, Phys. Rev. Lett. 101, 246810 (2008).
  * [14] E. Alba, X. F. -Gonzalvo, J. M. -Petit, J. K. Pachos, and J. J. G. -Ripoll, Phys. Rev. Lett. 107, 235301 (2011).
  * [15] L. Tarruell, D. Greif, T. Uehlinger, G. Jotzu and T. Esslinger, Nature 483, 302 (2012).
  * [16] G. Jotzu, M. Messer, R. Desbuquois, M. Lebrat, T. Uehlinger, D. Greif and T. Esslinger, Nature 515, 237–240 (2014).
  * [17] H. S. Kim and H. Y. Kee, npj Quant Mater 2, 20 (2017).
  * [18] Y. Ding, Y. Peng, Y. Zhu, X. Fan, J. Yang, B. Liang, X. Zhu, X. Wan, and J. Cheng, Phys. Rev. Lett. 122, 014302 (2019).
  * [19] K. F. Garrity and D. Vanderbilt, Phys. Rev. B 90, 121103(R) (2014).
  * [20] S. Spurrier and N. R. Cooper, Phys. Rev. Research 2, 033071 (2020).
  * [21] P. Cheng, P. W. Klein, K. Plekhanov, K. Sengstock, M. Aidelsburger, C. Weitenberg, and K. Le. Hur, Phys. Rev. B 100, 081107(R) (2019).
  * [22] S. Sorn, Phys. Rev. B 98, 125145 (2018).
  * [23] J. Panas, B. Irsigler, J.-H. Zheng, and W. Hofstetter, Phys. Rev. B 102, 075403 (2020).
  * [24] S. Mondal, P. Kapri, B. Dey, T. K. Ghosh and S. Basu, J. Phys.: Condens. Matter 33, 225504, (2021).
  * [25] S. Mondal and S. Basu, Phys. Rev. B 105, 235441 (2022).
  * [26] S. Mondal and S. Basu, Phys. Rev. B 107, 035421 (2023).
  * [27] V. Pardo and W.E. Pickett, Phys. Rev. Lett. 102, 166803 (2009).
  * [28] V. Pardo, W.E. Pickett, Phys. Rev. B 81, 035111 (2010).
  * [29] A. S. Rodin, A. Carvalho, and A. H. Castro Neto, Phys. Rev. Lett. 112, 176801 (2014).
  * [30] J. Guan, Z. Zhu, and D. Tománek, Phys. Rev. Lett. 113, 046804 (2014).
  * [31] Y. Suzumura, T. Morinari, and F. Piéchon,  J. Phys. Soc. Jpn. 82, 023708 (2013).
  * [32] Y. Hasegawa, R. Konno, H. Nakano, and M. Kohmoto, Phys. Rev. B 74, 033413 (2006).
  * [33] J. Kim, S. S. Baik, S. H. Ryu, Y. Sohn, S. Park, B.-G. Park, J. Denlinger, Y. Yi, H. J. Choi, and K. S. Kim, Science 349, 723 (2015).
  * [34] C. Zhong, Y. Chen, Y. Xie, Y.-Y. Sun, and S. Zhang, Phys. Chem. Chem. Phys. 19, 3820 (2017).
  * [35] D. Sticlet and F. Piéchon, Phys. Rev. B 87, 115402 (2013).
  * [36] S, Mondal and S. Basu, Phys. Rev. B 105, 235409 (2022).
  * [37] S. Yang, Z. -C. Gu, K. Sun, and S. D. Sarma, Phys. Rev. B 86, 241112(R) (2012).
  * [38] W. -C. Chen, R. Liu, Y. -F. Wang, and C. -D. Gong, Phys. Rev. B 86, 085311 (2012).
  * [39] Y. Yang, Y. F. Zhang, L. Sheng and D. Y. Xing, EPL 105, 27005 (2014).
  * [40] Y. Yang, X. Li, and D. Xing, Eur. Phys. J. B 89, 217 (2016).
  * [41] A. Alase and D. L. Feder, Phys. Rev. A 103, 053305 (2021).
  * [42] M. Lacki, J. Zarkrzewski and N. Goldman, SciPost Phys. 10, 112 (2021).
  * [43] C. Fang, M. J. Gilbert, and B. A. Bernevig, Phys. Rev. Lett. 112, 046801 (2014).
  * [44] J. Wang, B. Lian, H. Zhang, Y. Xu, and S.-C. Zhang Phys. Rev. Lett. 111, 136801 (2013).
  * [45] H. Zhao, T. Zhang, X. Zhang, M. -H. Lu, Y. -F. Chen,  Front. Phys. 10 (2022).
  * [46] W. Zhu, C. Song, H. Bai, L. Liao, and F. Pan Phys. Rev. B 105, 155122 (2022).
  * [47] J. Ge, Y. Liu, J. Li, H. Li, T. Luo, Y. Wu, Y. Xu, J. Wang National Science Review 8, 1280 (2020).
  * [48] Y. F. Zhao, R. Zhang, R. Mei, L. Zhou, H. Yi, Y. Q. Zhang, J. Yu, R. Xiao, K. Wang, N. Samarth, M. H. W. Chan, C. X. Liu, and C. Z. Chang, Nature (London) 588, 419 (2020).
  * [49] However, for obtaining the Chern number phase diagrams see Figs. 4, 5, and 8, $\phi_{l}$ and $\phi_{u}$ are varied in the range $[-\pi:+\pi]$.
  * [50] D. J. Thouless, Topological Quantum Numbers in Nonrelativistic Physics (World Scientific, Singapore, 1998).
  * [51] J. E. Avron, L. Sadun, J. Segert, and B. Simon, Phys. Rev. Lett. 61, 1329 (1988).
  * [52] C. -X. Liu, S. -C. Zhang, and X. -L. Qi Annu. Rev. Condens. Matter Phys. 7 301 (2016).
  * [53] D. Sticlet, F. Piéchon, J.-N. Fuchs, P. Kalugin, and P. Simon, Phys. Rev. B 85, 165456 (2012).
  * [54] K. Nakada, M. Fujita, G. Dresselhaus, and M. S. Dresselhaus, Phys. Rev. B 54, 17954 (1996).
  * [55] Y. Hatsugai Phys. Rev. B 48, 11851 (1988).
  * [56] D. Xiao, M.-C. Chang, and Q. Niu, Rev. Mod. Phys. 82, 1959 (2010).
  * [57] D. Culcer, A.MacDonald, and Q. Niu, Phys. Rev. B 68, 045327 (2003).