# Entropy- and Distance-Based Predictors From GPT-2 Attention Patterns Predict
Reading Times Over and Above GPT-2 Surprisal

Byung-Doh Oh
Department of Linguistics
The Ohio State University
<EMAIL_ADDRESS>
&William Schuler
Department of Linguistics
The Ohio State University
<EMAIL_ADDRESS>

###### Abstract

Transformer-based large language models are trained to make predictions about
the next word by aggregating representations of previous tokens through their
self-attention mechanism. In the field of cognitive modeling, such attention
patterns have recently been interpreted as embodying the process of cue-based
retrieval, in which attention over multiple targets is taken to generate
interference and latency during retrieval. Under this framework, this work
first defines an entropy-based predictor that quantifies the diffuseness of
self-attention, as well as distance-based predictors that capture the
incremental change in attention patterns across timesteps. Moreover, following
recent studies that question the informativeness of attention weights, we also
experiment with alternative methods for incorporating vector norms into
attention weights. Regression experiments using predictors calculated from the
GPT-2 language model show that these predictors deliver a substantially better
fit to held-out self-paced reading and eye-tracking data over a rigorous
baseline including GPT-2 surprisal. Additionally, the distance-based
predictors generally demonstrated higher predictive power, with effect sizes
of up to 6.59 ms per standard deviation on self-paced reading times (compared
to 2.82 ms for surprisal) and 1.05 ms per standard deviation on eye-gaze
durations (compared to 3.81 ms for surprisal).

##  1 Introduction

Much work in broad-coverage sentence processing has focused on studying the
role of expectation operationalized in the form of surprisal (Hale, 2001;
Levy, 2008) using language models (LMs) to define a conditional probability
distribution of a word given its context (Smith and Levy, 2013; Goodkind and
Bicknell, 2018). Recently, this has included Transformer-based language models
(Wilcox et al., 2020; Merkx and Frank, 2021; Oh et al., 2022).

Figure 1: The predictors and attention weight formulations examined in this
work. The entropy-based predictor (i.e. NAE) quantifies the diffuseness of
attention over previous tokens at a given timestep, while the distance-based
predictors (i.e. $\Delta$NAE, MD, EMD) capture the change in attention
patterns across consecutive timesteps (top row: weights at ‘journey,’ bottom
row: weights at ‘to’). These predictors can be calculated from attention
weights formulated using different methods (i.e. Attn-W, Attn-N, AttnRL-N).

However, expectation-based accounts have empirical shortcomings, such as being
unable to fully account for garden-path effects (van Schijndel and Linzen,
2021) or predict the timing of delays in certain constructions (Levy et al.,
2013). For this reason, some research has begun to focus on the effects of
memory and attention using predictors calculated from language model
representations. For example, Ryu and Lewis (2021) recently drew connections
between the self-attention patterns of Transformers (Vaswani et al., 2017) and
cue-based retrieval models of sentence comprehension (e.g. Lewis et al.,
2006). Their attention entropy, which quantifies the diffuseness of the
attention weights over previous tokens, showed patterns that are consistent
with similarity-based interference observed during the processing of subject-
verb agreement. However, these results relied on identifying one attention
head specialized for the nsubj dependency, and an aggregated version of this
predictor was not very strong in predicting naturalistic reading times in the
presence of a surprisal predictor (Ryu and Lewis, 2022).

This work therefore defines and evaluates several memory- and attention-based
predictors derived from the self-attention patterns of the Transformer-based
GPT-2 language model (Radford et al., 2019) on two naturalistic datasets, in
the presence of a strong GPT-2 surprisal baseline. First, normalized attention
entropy expands upon Ryu and Lewis’s (2021) attention entropy by re-
normalizing the attention weights and controlling for the number of tokens in
the previous context. Additionally, three distance-based predictors that
quantify the shift in attention patterns across consecutive timesteps are
presented, based on the idea that the reallocation of attentional focus
entails processing difficulty.

Moreover, motivated by work on interpreting large language models that
question the connection between attention weights and model predictions (e.g.
Jain and Wallace, 2019), the norm-based analysis of the transformed vectors
(Kobayashi et al., 2020, 2021) is newly applied to GPT-2 in this work to
inform alternative formulations of attention weights. For example, it has been
observed that while large language models tend to place high attention weights
on special tokens (e.g. <|endoftext|> of GPT-2), these tokens contribute very
little to final model predictions as their ‘value’ vectors have near-zero
norms (Kobayashi et al., 2020). Attention weight formulations that incorporate
the norms of the transformed vectors should therefore alleviate the over-
representation of such special tokens and represent the contribution of each
token more accurately.

Results from regression analyses using these predictors show significant and
substantial effects in predicting self-paced reading times and eye-gaze
durations during naturalistic reading, even in the presence of a robust
surprisal predictor. Additionally, alternative formulations of attention
weights that incorporate the norms of the transformed vectors are shown to
further improve the predictive power of these predictors.

##  2 Background

This section provides a mathematical definition of the self-attention
mechanism underlying the GPT-2 language model and describes alternative norm-
based formulations of attention weights.

###  2.1 Masked Self-Attention of GPT-2 Language Models

GPT-2 language models (Radford et al., 2019) use a variant of a multi-layer
Transformer decoder proposed in Vaswani et al. (2017). Each decoder layer
consists of a masked self-attention block and a feedforward neural network:

$\mathbf{x}_{l+1,i}=\text{\sc FF}(\text{\sc
LN}_{\text{out}}(\mathbf{o}_{l,i}+\mathbf{x}_{l,i}))+(\mathbf{o}_{l,i}+\mathbf{x}_{l,i}),$
(1)

where $\mathbf{x}_{l,i}\in\mathbb{R}^{d}$ is the $i$th input representation at
layer $l$, FF is a two-layer feedforward neural network, LN${}_{\text{out}}$
is a vector-wise layer normalization operation, and $\mathbf{o}_{l,i}$ is the
output representation from the multi-head self-attention mechanism, in which
$H$ multiple heads simultaneously mix representations from the previous
context. This output $\mathbf{o}_{l,i}$ can be decomposed into the sum of
representations resulting from each attention head $h$:

$\mathbf{o}_{l,i}=\sum_{h=1}^{H}{\mathbf{V}_{h}\begin{bmatrix}\mathbf{X^{\prime}}_{l,i}\\\
\mathbf{1}^{\top}\end{bmatrix}\mathbf{a}_{l,h,i}},$ (2)

where $\mathbf{X^{\prime}}_{l,i}\overset{\rm
def}{=}[\mathbf{x^{\prime}}_{l,1},...,\mathbf{x^{\prime}}_{l,i}]\in\mathbb{R}^{d\times
i}$ is the sequence of layer-normalized input representations leading up to
$\mathbf{x^{\prime}}_{l,i}$ from the previous context, and
$\mathbf{x^{\prime}}_{l,i}\overset{\rm def}{=}\text{\sc
LN}_{\text{in}}(\mathbf{x}_{l,i})$ is the layer-normalized version of
$\mathbf{x}_{l,i}$. $\mathbf{V}_{h}$ represents the head-specific value-output
transformation,111As $\begin{bmatrix}\mathbf{W}\>\mathbf{b}\end{bmatrix}$
$\begin{bmatrix}\mathbf{x}\\\ 1\end{bmatrix}=\mathbf{W}\mathbf{x}+\mathbf{b}$,
the bias vectors are omitted from the equations. Additionally, for the
simplicity of notation, the ‘output’ affine transform, which applies to the
concatenated ‘value-transformed’ vectors from all attention heads in a typical
implementation, is subsumed into $\mathbf{V}_{h}$. The bias vector of the
‘output’ transform is assumed to be distributed equally across heads. and
$\mathbf{a}_{l,h,i}\in\mathbb{R}^{i}$ is the vector of attention weights:

$\mathbf{a}_{l,h,i}=\underset{}{\text{\sc
SoftMax}}\bigg{(}\frac{\big{(}\mathbf{K}_{h}\\!\begin{bmatrix}\mathbf{X^{\prime}}_{l,i}\\\
\mathbf{1}^{\top}\end{bmatrix}\big{)}^{\top}\mathbf{Q}_{h}\\!\begin{bmatrix}\mathbf{x^{\prime}}_{l,i}\\\
1\end{bmatrix}}{\sqrt{d_{h}}}\bigg{)},$ (3)

where $\mathbf{Q}_{h}$ and $\mathbf{K}_{h}$ represent the head-specific query
and key transformations respectively, and $d_{h}=d/H$ is the dimension of each
attention head.

Figure 2: Computations performed within the self-attention block of one head
of the GPT-2 language model at a given timestep ($i\\!=\\!5$). While the
masked self-attention mechanism aggregates representations from the previous
context in a typical implementation, the linear nature of the subsequent
computations allows this aggregation to be deferred to after the residual
connection and layer normalization, thereby allowing updated representations
to inform alternative formulations of attention weights (i.e. Attn-N,
AttnRL-N).

LNα is a vector-wise layer normalization operation (Ba et al., 2016) that
first standardizes the vector and subsequently conducts elementwise
transformations using learnable parameters
$\mathbf{c}_{\alpha},\mathbf{b}_{\alpha}\in\mathbb{R}^{d}$:

$\text{\sc
LN}_{\alpha}(\mathbf{y})=\frac{\mathbf{y}-m(\mathbf{y})}{s(\mathbf{y})}\odot\mathbf{c}_{\alpha}+\mathbf{b}_{\alpha},$
(4)

where $\alpha\in\\{\text{in},\text{out}\\}$, $m(\mathbf{y})$ and
$s(\mathbf{y})$ denote the elementwise mean and standard deviation
respectively, and $\odot$ denotes a Hadamard product.

###  2.2 Weight- and Norm-Based Analysis of Previous Context

Previous work that has studied the inner mechanism of Transformer models has
focused on analyzing the relative contribution of each token to its final
prediction. As a measure that quantifies the ‘strength’ of contribution,
attention weights from the self-attention mechanism have been most commonly
used. Similarly to recent work in cognitive modeling (e.g. Ryu and Lewis,
2021), this work also evaluates predictors calculated from attention weights
(Attn-W).

$\mathbf{a}_{\text{W},l,h,i}=\mathbf{a}_{l,h,i}$ (5)

While analysis using attention weights is common, the assumption that
attention weights alone reflect the contribution of each token disregards the
magnitudes of the transformed input vectors, as pointed out by Kobayashi et
al. (2020). As an alternative, they proposed a norm-based analysis of the
self-attention mechanism, which quantifies the contribution of each token as
the norm of the transformed vector multiplied by the attention weight. In this
work, in order to quantify the relative contribution of each token in the
previous context, the norms of the transformed vectors are normalized across
the sequence, resulting in ‘norm-aware’ weights that are comparable to
attention weights (Attn-N).

$\displaystyle\mathbf{a}_{\text{N},l,h,i}=$
$\displaystyle\medmath{\frac{\bigg{(}\bigg{(}\mathbf{V}_{h}\\!\begin{bmatrix}\mathbf{X^{\prime}}_{l,i}\\\
\mathbf{1}^{\top}\end{bmatrix}\odot\mathbf{V}_{h}\\!\begin{bmatrix}\mathbf{X^{\prime}}_{l,i}\\\
\mathbf{1}^{\top}\end{bmatrix}\bigg{)}^{\\!\\!\top}\mathbf{1}\bigg{)}^{\\!\frac{1}{2}}\\!\odot\mathbf{a}_{l,h,i}}{\mathbf{1}^{\top}\bigg{(}\\!\bigg{(}\\!\bigg{(}\mathbf{V}_{h}\\!\begin{bmatrix}\mathbf{X^{\prime}}_{l,i}\\\
\mathbf{1}^{\top}\end{bmatrix}\odot\mathbf{V}_{h}\\!\begin{bmatrix}\mathbf{X^{\prime}}_{l,i}\\\
\mathbf{1}^{\top}\end{bmatrix}\bigg{)}^{\\!\\!\top}\mathbf{1}\bigg{)}^{\\!\frac{1}{2}}\\!\odot\mathbf{a}_{l,h,i}\bigg{)}}}$
(6)

More recently, Kobayashi et al. (2021) showed that the residual connection and
the layer normalization operation (ResLN; RL) that follow the self-attention
mechanism can also be conducted before aggregating representations over token
positions. Motivated by this observation, the vector norms that take into
consideration these subsequent operations are also examined in this work.
Similarly to Attn-N, the norms are normalized across the sequence to yield
weights that are comparable (AttnRL-N):

$\displaystyle\mathbf{a}_{\text{RL-N},l,h,i}=$
$\displaystyle\medmath{\frac{((g(\mathbf{V}_{h}\\!\begin{bmatrix}\mathbf{X^{\prime}}_{l,i}\\\
\mathbf{1}^{\top}\end{bmatrix})\odot
g(\mathbf{V}_{h}\\!\begin{bmatrix}\mathbf{X^{\prime}}_{l,i}\\\
\mathbf{1}^{\top}\end{bmatrix}))^{\top}\mathbf{1})^{\frac{1}{2}}\odot\mathbf{a}_{l,h,i}}{\mathbf{1}^{\top}(((g(\mathbf{V}_{h}\\!\begin{bmatrix}\mathbf{X^{\prime}}_{l,i}\\\
\mathbf{1}^{\top}\end{bmatrix})\odot
g(\mathbf{V}_{h}\\!\begin{bmatrix}\mathbf{X^{\prime}}_{l,i}\\\
\mathbf{1}^{\top}\end{bmatrix}))^{\top}\mathbf{1})^{\frac{1}{2}}\\!\odot\\!\mathbf{a}_{l,h,i})}},$
(7)

where $g(\cdot)$ incorporates the residual connection ($+\mathbf{x}_{l,i}$)
and layer normalization ($\text{LN}_{\text{out}}$) of Eq. 1. Following the
assumption that the residual connection serves to ‘preserve’ the
representation at position $i$ (Kobayashi et al., 2021) and that it is
distributed equally across heads, $\mathbf{x}_{l,i}$ is added to the
representation at position $i$ of each head after dividing it by the number of
heads:

$\displaystyle g(\mathbf{Y})\overset{\rm def}{=}$
$\displaystyle\medmath{\sum_{j=1}^{i}\\!\left(\\!\\!\begin{cases}\\!\frac{\mathbf{Y}\delta_{j}+\frac{\mathbf{x}_{l,i}}{H\mathbf{a}_{l\\!,\\!h\\!,\\!i\\![\\!i\\!]}}-m(\mathbf{Y}\delta_{j}+\frac{\mathbf{x}_{l,i}}{H\mathbf{a}_{l\\!,\\!h\\!,\\!i\\![\\!i\\!]}}\\!)}{s(\mathbf{o}_{l,i}+\mathbf{x}_{l,i})}\\!\odot\\!\mathbf{c}_{\text{out}}\\!+\\!\frac{\mathbf{b}_{\text{out}}}{H}&\\!\\!\\!\\!\text{if
}i\\!=\\!j\\\
\\!\frac{\mathbf{Y}\delta_{j}-m(\mathbf{Y}\delta_{j})}{s(\mathbf{o}_{l,i}+\mathbf{x}_{l,i})}\\!\odot\\!\mathbf{c}_{\text{out}}\\!+\\!\frac{\mathbf{b}_{\text{out}}}{H}&\\!\\!\\!\\!\text{if
}i\\!\neq\\!j\end{cases}\\!\\!\right)\\!\\!\delta_{j}^{\top}}\\!,$ (8)

where $\delta_{j}$ is a Kronecker delta vector consisting of a one at element
$j$ and zeros elsewhere, and $\mathbf{c}_{\text{out}}$ and
$\mathbf{b}_{\text{out}}$ refer to the learnable parameters of $\text{\sc
LN}_{\text{out}}$.

##  3 Entropy- and Distance-Based Predictors From Attention Patterns

Given the different formulations of self-attention weights, entropy-based
predictors that quantify the diffuseness of self-attention and distance-based
predictors that capture the incremental change in attention patterns across
timesteps can be defined. The first predictor defined in this work is
normalized attention entropy (NAE):

$\displaystyle\text{\sc NAE}_{\pi,l,h,i}=$
$\displaystyle\medmath{\frac{\mathbf{a}_{\pi,l,h,i[1:i-1]}^{\top}}{\log_{2}\\!{(i-1)}\mathbf{1}^{\top}\mathbf{a}_{\pi,l,h,i[1:i-1]}}(\log_{2}\\!{\frac{\mathbf{a}_{\pi,l,h,i[1:i-1]}}{\mathbf{1}^{\top}\mathbf{a}_{\pi,l,h,i[1:i-1]}}})},$
(9)

where $\pi\in\\{\text{W},\text{N},\text{RL-N}\\}$. This is similar to the
attention entropy proposed by Ryu and Lewis (2021) as a measure of
interference in cue-based recall attributable to uncertainty about the target,
with two notable differences. First, NAE controls for the number of tokens in
the previous context by normalizing the entropy by the maximum entropy that
can be achieved at timestep $i$. Furthermore, NAE also uses weights over
$\mathbf{x^{\prime}}_{l,1}$, …, $\mathbf{x^{\prime}}_{l,i-1}$ that have been
re-normalized to sum to 1, thereby adhering closer to the definition of
entropy, in which the mass of interest sums to 1.222Preliminary analyses
showed that regression models with attention entropy proper without these
adjustments failed to converge coherently.

In addition to NAE, distance-based predictors are defined for capturing
effortful change in attention patterns across timesteps. However, as it
currently remains theoretically unclear how this distance should be defined,
this exploratory work sought to provide empirical results for different
distance functions. The first is $\Delta$NAE, which quantifies the change in
diffuseness across timesteps:

$\Delta\text{\sc NAE}_{\pi,l,h,i}=|\text{\sc NAE}_{\pi,l,h,i}\\!-\\!\text{\sc
NAE}_{\pi,l,h,i-1}|$ (10)

As with NAE, this predictor is insensitive to how the attention weights are
reallocated between tokens in the previous context to the extent that the
overall diffuseness remains unchanged.

The second distance-based predictor is Manhattan distance (MD).333For the
purpose of calculating this predictor, the $i$th element of
$\mathbf{a}_{\pi,l,h,i-1}$ is assumed to be 0.

$\text{\sc
MD}_{\pi,l,h,i}=||\mathbf{a}_{\pi,l,h,i}-\mathbf{a}_{\pi,l,h,i-1}||_{1}$ (11)

MD directly measures the magnitude of change in attention weights over all
tokens at timestep $i$. MD is less sensitive to the linear distance between
tokens and therefore makes it consistent with the predictions of McElree et
al. (2003), who found that processing speed was not influenced by the amount
of intervening linguistic material in the formation of a dependency.

Finally, the Earth Mover’s Distance (EMD; Rubner et al., 2000) is applied to
quantify the shift in attention weights. EMD is derived from a solution to the
Monge-Kantorovich problem (Rachev, 1985), which aims to minimize the amount of
“work” necessary to transform one histogram into another. Formally, let
$P=\\{(p_{1},w_{p_{1}}),...,(p_{m},w_{p_{m}})\\}$ be the first histogram with
$m$ bins, where $p_{r}$ represents the bin and $w_{p_{r}}$ represents the
weight of the bin; $Q=\\{(q_{1},w_{q_{1}}),...,(q_{n},w_{q_{n}})\\}$ the
second histogram with $n$ bins; and $\mathbf{D}=[d_{rs}]$ the distance matrix
where $d_{rs}$ is the ground distance between bins $p_{r}$ and $q_{s}$. The
problem is to find an optimal flow $\mathbf{F}=[f_{rs}]$, where $f_{rs}$
represents the flow between $p_{r}$ and $q_{s}$, that minimizes the overall
work.

$\text{\sc Work}(P,Q,\mathbf{F})=\sum_{r=1}^{m}\sum_{s=1}^{n}{d_{rs}f_{rs}}$
(12)

Once the optimal flow is found, the EMD is defined as the work normalized by
the total flow.444The optimal flow can be found using the transportation
simplex method. Additionally, due to the constraint that the total flow is
equal to $\text{min}(\sum_{r=1}^{m}{w_{p_{r}}},\sum_{s=1}^{n}{w_{q_{s}}}$),
the total flow is always equal to 1 in the context of attention weights.

$\text{\sc
EMD}(P,Q)=\frac{\sum_{r=1}^{m}\sum_{s=1}^{n}d_{rs}f_{rs}}{\sum_{r=1}^{m}\sum_{s=1}^{n}f_{rs}}$
(13)

To quantify the minimum amount of work necessary to ‘transform’ the attention
weights, the EMD between attention weights at consecutive timesteps is
calculated using EMD($\mathbf{a}_{\pi,l,h,i-1}$, $\mathbf{a}_{\pi,l,h,i}$).
The ground distance is defined as $d_{rs}=\frac{|r-s|}{i-1}$ in order to
control for the number of tokens in the previous context. EMD can be
interpreted as being consistent with Dependency Locality Theory (Gibson, 2000)
in that reallocating attention weights to tokens further away in the context
incurs more cost than reallocating weights to closer tokens.

Code for calculating all of the predictors from GPT-2 under the different
attention weight formulations is publicly available at
https://github.com/byungdoh/attn_dist.

##  4 Experiment 1: Evaluation of Predictors on Human Reading Times

In order to evaluate the contribution of the entropy- and distance-based
predictors, regression models containing commonly used baseline predictors,
surprisal predictors, and one predictor of interest were fitted to self-paced
reading times and eye-gaze durations collected during naturalistic language
processing. In this work, we adopt a statistical procedure that directly
models temporal diffusion (i.e. a lingering reponse to stimuli) by estimating
continuous impulse response functions and controls for overfitting by
assessing the external validity of these predictors through a non-parametric
test on held-out data.

###  4.1 Response Data

The first experiment described in this paper used the Natural Stories Corpus
(Futrell et al., 2021), which contains self-paced reading times from 181
subjects that read 10 naturalistic stories consisting of 10,245 words. The
data were filtered to exclude observations for sentence-initial and sentence-
final words, observations from subjects who answered fewer than four
comprehension questions correctly, and observations with durations shorter
than 100 ms or longer than 3000 ms. This resulted in a total of 770,102
observations, which were subsequently partitioned into a fit partition of
384,905 observations, an exploratory partition of 192,772 observations, and a
held-out partition of 192,425 observations.555For both datasets, the fit
partition, exploratory partition, and held-out partition contain data points
whose summed subject and sentence number have modulo four equal to zero or
one, two, and three respectively. The partitioning allows model selection
(e.g. making decisions about baseline predictors and random effects structure)
to be conducted on the exploratory partition and a single hypothesis test to
be conducted on the held-out partition, thus obviating the need for multiple
trials correction. All observations were log-transformed prior to regression
modeling.

Additionally, the set of go-past durations from the Dundee Corpus (Kennedy et
al., 2003) also provided the response variable for regression modeling. The
Dundee Corpus contains eye-gaze durations from 10 subjects that read 67
newspaper editorials consisting of 51,501 words. The data were filtered to
remove unfixated words, words following saccades longer than four words, and
words at sentence-, screen-, document-, and line-starts and ends. This
resulted in a total of 195,507 observations, which were subsequently
partitioned into a fit partition of 98,115 observations, an exploratory
partition of 48,598 observations, and a held-out partition of 48,794
observations. All observations were log-transformed before model fitting.

###  4.2 Predictors

For each dataset, a set of baseline predictors that capture basic, low-level
cognitive processing were included in all regression models.

  * •

Self-paced reading times (Futrell et al., 2021): word length measured in
characters, index of word position within each sentence;

  * •

Eye-gaze durations (Kennedy et al., 2003): word length measured in characters,
index of word position within each sentence, saccade length, whether or not
the previous word was fixated.

In addition to the baseline predictors, two surprisal predictors were also
included in all regression models evaluated in this experiment. The first is
unigram surprisal as a measure of word frequency, which was calculated using
the KenLM toolkit (Heafield et al., 2013) with parameters estimated on the
English Gigaword Corpus (Parker et al., 2009). The second is surprisal from
GPT-2 Small (Radford et al., 2019), which is trained on $\sim$8B tokens of the
WebText dataset. Surprisal from the smallest GPT-2 model was chosen because it
has been shown to be the most predictive of self-paced reading times and eye-
gaze durations among surprisal from all variants of GPT-2 (Oh et al., 2022).

Finally, the entropy- and distance-based predictors defined in Section 3 were
calculated from the attention patterns (i.e. $\mathbf{a}_{\pi,l,h,i}$ where
$\pi~{}\in~{}\\{\text{W},\text{N},\text{RL-N}\\}$) of heads on the topmost
layer of GPT-2 Small. Contrary to previous studies that analyzed the attention
patterns of all layers, this work focuses on analyzing those of the topmost
layer out of the concern that the attention patterns of lower layers are less
interpretable to the extent that they perform intermediate computations for
the upper layers. Since the topmost layer generates the representation that is
used for model prediction, the attention patterns from this layer are assumed
to reflect the contribution of each previous token most directly.
Subsequently, the by-head predictors were aggregated across heads to calculate
by-word predictors. This assumes that each attention head contributes equally
to model prediction, and is also consistent with the formulation of multi-head
self-attention in Eq. 2.

Figure 3: Improvements in CDR model log-likelihood from including each
predictor on the exploratory partition of Natural Stories self-paced reading
data (left) and Dundee eye-tracking data (right).

To calculate surprisal as well as the entropy- and distance-based predictors,
each story of the Natural Stories Corpus and each article of the Dundee Corpus
was tokenized according GPT-2’s byte-pair encoding (BPE; Sennrich et al.,
2016) tokenizer and was provided as input to the GPT-2 model. In cases where
each story/article did not fit into a single context window, the second half
of the previous context window served as the first half of a new context
window to calculate predictors for the remaining tokens.666In practice, most
stories/articles fit within two windows. Additionally, when a single word
$w_{i}$ was tokenized into multiple subword tokens, negative log probabilities
of subword tokens corresponding to $w_{i}$ were added together to calculate
$\mathsf{S}(w_{i})=-\log\mathsf{P}(w_{i}\mid w_{1..i-1})$. Similarly, the
entropy- and distance-based predictors for such subword tokens were also added
together.

###  4.3 Procedures

To evaluate the predictive power of each predictor of interest, a baseline
regression model containing only the baseline predictors and full regression
models containing one predictor of interest on top of the baseline regression
model were first fitted to the fit partition of each dataset. In order to
control for the confound of temporal diffusion, continuous impulse response
functions (IRFs) were estimated using the statistical framework of continuous-
time deconvolutional regression (CDR; Shain and Schuler, 2021). The predictors
of interest were centered, and all regression models included a by-subject
random intercept.777Early analysis on the exploratory partition revealed that
the training data does not support a richer random-effects structure, which
led to severe overfitting. Please refer to Appendix A for details on IRF
specifications, the optimization procedure of CDR models, as well as the
transformations of baseline predictors.

As a preliminary analysis, the predictive power of different predictors was
compared on the exploratory partition by calculating the increase in log-
likelihood ($\Delta$LL) to the baseline regression model as a result of
including the predictor, following recent work (Goodkind and Bicknell, 2018;
Wilcox et al., 2020; Oh et al., 2021). Subsequently, based on the preliminary
exploratory results, the predictive power of one best entropy-based predictor
and that of one best distance-based predictor were evaluated on the held-out
partition of both datasets. More specifically, statistical significance
testing was conducted using a paired permutation test (Demšar, 2006) of the
difference in by-item squared error between the baseline regression model and
the respective full regression models.

###  4.4 Results

The results in Figure 3 show that across both corpora, most of the entropy-
and distance-based predictors make a notable contribution to regression model
fit under all attention formulations. Given that the baseline model contains
strong predictors such as unigram surprisal and GPT-2 surprisal, this may
suggest their validity as predictors of comprehension difficulty. The
exception to this is the EMD predictor, which did not show an increase in
likelihood under Attn-W on the Natural Stories Corpus and under AttnRL-N on
the Dundee Corpus. As EMD is sensitive to how the ground distance $d_{rs}$
between bins is defined (i.e. how costly it is to reallocate attention weights
across input positions), a more principled definition of $d_{rs}$ may make EMD
a more robust predictor.888For example, a recency bias may be incorporated by
defining the ground distance from tokens closer to the current timestep to be
smaller than that from tokens that are farther back in the previous context.
Across the two corpora, Attn-N+NAE and AttnRL-N+MD appear to be the most
predictive among the entropy- and distance-based predictors respectively.

Additionally, the NAE and $\Delta$NAE predictors showed different trends
across the two corpora, where NAE contributed to stronger model fit than
$\Delta$NAE on the Natural Stories Corpus, while the opposite trend was
observed on the Dundee Corpus. In contrast to various surprisal predictors
that have shown a very similar trend in terms of predictive power across these
two corpora (Oh et al., 2022), these two predictors may shed light on the
subtle differences between self-paced reading times and eye-gaze durations.
Finally, incorporating vector norms into attention weights (i.e. Attn-N and
AttnRL-N) generally seems to improve the predictive power of these predictors,
which provides support for the informativeness of input vectors in analyzing
attention patterns (Kobayashi et al., 2020, 2021).

Table 1 presents the effect sizes of Attn-N+NAE and AttnRL-N+MD on the held-
out partition of the Natural Stories Corpus and the Dundee Corpus, which were
derived by calculating how much increase in reading times the regression model
would predict at average predictor value given an increase of one standard
deviation. On both datasets, AttnRL-N+MD appears to be a strong predictor of
reading times, which contributed to significantly lower held-out errors. The
entropy-based Attn-N+NAE predictor showed contrasting results across the two
corpora, showing a large effect size on the Natural Stories Corpus but not on
the Dundee Corpus. This is consistent with the differential results between
NAE and $\Delta$NAE on the exploratory partition of the two corpora and may
hint at differences between self-paced reading times and eye-gaze durations.
In terms of magnitude, the two predictors showed large effect sizes on the
Natural Stories Corpus, which were more than twice that of GPT-2 surprisal. On
the Dundee Corpus, however, the effect size of AttnRL-N+MD was much smaller
compared to that of GPT-2 surprisal.

Corpus | Predictor | Effect Size ($p$-value)
---|---|---
Natural Stories | Attn-N+NAE | 6.87 ms ($p<0.001$)
GPT2Surp | 2.56 ms
| AttnRL-N+MD | 6.59 ms ($p<0.001$)
| GPT2Surp | 2.82 ms
Dundee | Attn-N+NAE | N/A999The estimated IRF for Attn-N+NAE showed a sign error, likely due to poor convergence. Therefore, we treated the effect of this predictor to be statistically non-significant. (n.s.)
GPT2Surp | 4.22 ms
| AttnRL-N+MD | 1.05 ms ($p<0.001$)
| GPT2Surp | 3.81 ms
Table 1: The per standard deviation effect sizes of the predictors on the
held-out partition of the Natural Stories Corpus and the Dundee Corpus.
Statistical significance was determined by a paired permutation test of the
difference in by-item squared error between the baseline regression model and
the respective full regression model containing the predictor of interest. The
effect sizes of GPT-2 surprisal from the same regression models are presented
for comparison.

##  5 Experiment 2: Do NAE and MD Independently Explain Reading Times?

The previous experiment revealed that on the Natural Stories Corpus, the
select entropy- and distance-based predictors from the attention patterns of
GPT-2 contributed to significantly higher regression model fit. Although they
showed similarly large effect sizes on the held-out partition, the two
predictors may independently explain reading times, as they are defined to
quantify different aspects of attention patterns. The second experiment
examines this possibility following similar procedures as the previous
experiment.

###  5.1 Procedures

In order to determine whether the effect of one predictor subsumes that of the
other, a CDR model including both Attn-N+NAE and AttnRL-N+MD was first fit to
self-paced reading times of the fit partition of the Natural Stories Corpus.
The CDR model followed the same specifications, random effects structure, and
baseline predictors as described in Experiment 1. Subsequently, the fit of
this regression model on the held-out partition was compared to those of the
two regression models that contain only one of the two predictors from the
previous experiment. More specifically, the $\Delta$LL as a result of
including the predictor(s) of interest were calculated for each regression
model, and statistical significance testing was conducted using a paired
permutation test of the difference in by-item squared error between the new
‘Attn-N+NAE & AttnRL-N+MD’ regression model and the respective ‘Attn-N+NAE’
and ‘AttnRL-N+MD’ regression models, which allowed the contribution of each
predictor to be analyzed.

Figure 4: Improvements in CDR model log-likelihood from including Attn-N+NAE,
AttnRL-N+MD, and both on the held-out partition of Natural Stories self-paced
reading data.

###  5.2 Results

The results in Figure 4 show that the regression model including both
Attn-N+NAE and AttnRL-N+MD achieves higher $\Delta$LL on the held-out
partition of the Natural Stories Corpus compared to regression models
including only one of these predictors. Moreover, the difference in by-item
squared error between the ‘Attn-N+NAE & AttnRL-N+MD’ regression model and the
‘AttnRL-N+MD’ regression model was statistically significant ($p<0.001$). In
contrast, the significance between the ‘Attn-N+NAE & AttnRL-N+MD’ regression
model and the ‘Attn-N+NAE’ regression model was not statistically significant
($p=0.236$). This indicates that the entropy-based Attn-N+NAE contributes to
model fit over and above the distance-based AttnRL-N+MD and also subsumes its
effect in predicting reading times.

Figure 5: Pearson correlation coefficients between the predictors examined in
this work on the Natural Stories Corpus (upper right triangle) and the Dundee
Corpus (lower left triangle).

##  6 Correlation Analysis Between Predictors and Syntactic Categories

###  6.1 Procedures

In order to shed more light on the predictors newly proposed in this work, a
series of correlation analyses were conducted. First, Pearson correlation
coefficients were calculated between the entropy- and distance-based
predictors as well as the surprisal predictors to examine the similarity
between predictors and the influence of different attention weight
formulations. Subsequently, the most predictive Attn-N+NAE and AttnRL-N+MD
predictors were analyzed with a focus on identifying potential linguistic
correlates. This analysis used a version of the Natural Stories Corpus and the
Dundee Corpus that had been manually annotated using a generalized categorial
grammar annotation scheme (Shain et al., 2018).

###  6.2 Results

The correlation matrix in Figure 5 shows that within the same predictor, the
different attention formulations did not make a very large difference, with
‘intra-predictor’ correlation coefficients at above $0.85$ for most
predictors. An exception to this trend was EMD, which showed a correlation
coefficient of $0.74$ on the Natural Stories Corpus and $0.72$ on the Dundee
Corpus between Attn-W+EMD and AttnRL-N+EMD. Such difference is to be expected,
as EMD is the most sensitive to the change of location in attention weights
among the distance-based predictors. This is also consistent with the
exploratory regression results in Figure 3, where the $\Delta$LL of EMD
predictors varied the most as a function of different attention weight
formulations.

Additionally, the norm-based attention formulations seem to bring the entropy-
based NAE closer to the distance-based predictors in terms of correlation
coefficients. On both corpora, Attn-N+NAE and AttnRL-N+NAE show stronger
correlations to distance-based predictors than Attn-W+NAE. Interestingly, the
highest correlation coefficient between NAE and a distance-based predictor is
observed between Attn-N+NAE and AttnRL-N+MD at $0.90$ on the Natural Stories
Corpus and $0.91$ on the Dundee Corpus, which were the two strongest
predictors identified in Experiment 1. Such high correlation partially
explains the results of Experiment 2, in which the influence of Attn-N+NAE
subsumed that of AttnRL-N+MD.

Finally, the entropy- and distance-based predictors showed moderate
correlation to unigram surprisal at around $0.5$ on both corpora. With regard
to GPT-2 surprisal, these predictors showed weak correlation at around $0.3$
on the Natural Stories Corpus, and around $0.4$ on the Dundee Corpus. Together
with the regression results, this further suggests that the proposed
predictors capture a mechanistic process that is distinct from the frequency
or predictability of the word.

An analysis of the predictors according to syntactic categories showed that
AttnRL-N+MD may be sensitive to the transition between the subject constituent
and the predicate constituent. Its average values for nouns in different
contexts presented in Table 2 show that on both corpora, nouns occurring at
the end of subject NPs were associated with a greater shift in attention
patterns. This is consistent with the intuition that transitions between
constituents entails cognitive effort during incremental processing.

Corpus | Nouns | Count | Mean
---|---|---|---
Natural Stories | End of subj. NP | 1,051 | 19.024
Other | 3,098 | 16.520
Dundee | End of subj. NP | 4,496 | 18.982
Other | 15,763 | 16.661
Table 2: Average AttnRL-N+MD values for nouns at the end of subject NPs and
nouns in other contexts.

##  7 Discussion and Conclusion

This work builds on recent efforts to derive memory- and attention-based
predictors of processing difficulty to complement expectation-based accounts
of sentence processing. Based on the observation that the self-attention
patterns of Transformer-based language models can be interpreted as embodying
the process of cue-based retrieval, an entropy-based predictor that quantifies
the diffuseness of self-attention was first defined. Moreover, based on the
idea that reallocation of attention may incur processing difficulty, distance-
based predictors that capture the incremental change in attention patterns
across timesteps were defined. Regression results using these predictors
calculated from the GPT-2 language model showed that these entropy- and
distance-based predictors deliver a substantially better fit to self-paced
reading and eye-tracking data over a strong baseline including GPT-2
surprisal.

To our knowledge, this work is the first to report robust effects of
Transformer attention-based predictors in predicting reading times of broad-
coverage naturalistic data. This provides support for Ryu and Lewis’s (2021)
observation that the self-attention mechanism of Transformers embodies the
process of cue-based retrieval, and further suggests that representations that
exhibit similarity-based interference can be learned from the self-supervised
next-word prediction task. Additionally, the strength of the distance-based
predictors further demonstrates the potential to bring together expectation-
and memory-based theories of sentence processing under a coherent framework.

## Acknowledgments

We thank the reviewers for their helpful comments. This work was supported by
the National Science Foundation grant #1816891. All views expressed are those
of the authors and do not necessarily reflect the views of the National
Science Foundation.

## Limitations

The connection between attention patterns of Transformer-based language models
and human sentence processing drawn in this work is based on a model trained
on English text and data from human subjects that are native speakers of
English. Therefore, the connection made in this work may not generalize to
other languages. Additionally, although the alternative formulations of self-
attention weights resulted in stronger predictors of processing difficulty,
they are more computationally expensive to calculate as they rely on an
explicit decomposition of the matrix multiplication operation, which are
highly optimized in most packages.

## Ethics Statement

Experiments presented in this work used datasets from previously published
research (Futrell et al., 2021; Kennedy et al., 2003), in which the procedures
for data collection and validation are outlined. As this work focuses studying
the possible connection between attention patterns of large language models
and human sentence processing, its potential negative impacts on society seem
to be minimal.

## References

  * Ba et al. (2016) Jimmy Lei Ba, Jamie Ryan Kiros, and Geoffrey E. Hinton. 2016.  Layer normalization.  _arXiv preprint_.
  * Demšar (2006) Janez Demšar. 2006.  Statistical comparisons of classifiers over multiple data sets.  _Journal of Machine Learning Research_ , 7(1):1–30.
  * Dozat (2016) Timothy Dozat. 2016.  Incorporating Nesterov momentum into Adam.  In _ICLR Workshop Track_.
  * Futrell et al. (2021) Richard Futrell, Edward Gibson, Harry J. Tily, Idan Blank, Anastasia Vishnevetsky, Steven Piantadosi, and Evelina Fedorenko. 2021.  The Natural Stories corpus: A reading-time corpus of English texts containing rare syntactic constructions.  _Language Resources and Evaluation_ , 55:63–77.
  * Gibson (2000) Edward Gibson. 2000.  The Dependency Locality Theory: A distance-based theory of linguistic complexity.  In _Image, language, brain: Papers from the first mind articulation project symposium_ , pages 95–126, Cambridge, MA. MIT Press.
  * Goodkind and Bicknell (2018) Adam Goodkind and Klinton Bicknell. 2018.  Predictive power of word surprisal for reading times is a linear function of language model quality.  In _Proceedings of the 8th Workshop on Cognitive Modeling and Computational Linguistics_ , pages 10–18.
  * Hale (2001) John Hale. 2001.  A probabilistic Earley parser as a psycholinguistic model.  In _Proceedings of the Second Meeting of the North American Chapter of the Association for Computational Linguistics on Language Technologies_ , pages 1–8.
  * Heafield et al. (2013) Kenneth Heafield, Ivan Pouzyrevsky, Jonathan H. Clark, and Philipp Koehn. 2013.  Scalable modified Kneser-Ney language model estimation.  In _Proceedings of the 51st Annual Meeting of the Association for Computational Linguistics_ , pages 690–696.
  * Jain and Wallace (2019) Sarthak Jain and Byron C. Wallace. 2019.  Attention is not explanation.  In _Proceedings of the 2019 Conference of the North American Chapter of the Association for Computational Linguistics: Human Language Technologies, Volume 1 (Long and Short Papers)_ , pages 3543–3556.
  * Kennedy et al. (2003) Alan Kennedy, Robin Hill, and Joël Pynte. 2003.  The Dundee Corpus.  In _Proceedings of the 12th European Conference on Eye Movement_.
  * Kingma and Ba (2015) Diederik P. Kingma and Jimmy Ba. 2015.  Adam: A method for stochastic optimization.  In _Proceedings of the 3rd International Conference on Learning Representations_.
  * Kobayashi et al. (2020) Goro Kobayashi, Tatsuki Kuribayashi, Sho Yokoi, and Kentaro Inui. 2020.  Attention is not only a weight: Analyzing transformers with vector norms.  In _Proceedings of the 2020 Conference on Empirical Methods in Natural Language Processing_ , pages 7057–7075.
  * Kobayashi et al. (2021) Goro Kobayashi, Tatsuki Kuribayashi, Sho Yokoi, and Kentaro Inui. 2021.  Incorporating residual and normalization layers into analysis of masked language models.  In _Proceedings of the 2021 Conference on Empirical Methods in Natural Language Processing_ , pages 4547–4568.
  * Levy (2008) Roger Levy. 2008.  Expectation-based syntactic comprehension.  _Cognition_ , 106(3):1126–1177.
  * Levy et al. (2013) Roger Levy, Evalina Fedorenko, and Edward Gibson. 2013.  The syntactic complexity of Russian relative clauses.  _Journal of Memory and Language_ , 69(4):461–495.
  * Lewis et al. (2006) Richard L. Lewis, Shravan Vasishth, and Julie A. Van Dyke. 2006.  Computational principles of working memory in sentence comprehension.  _Trends in Cognitive Science_ , 10(10):447–454.
  * McElree et al. (2003) Brian McElree, Stephani Foraker, and Lisbeth Dyer. 2003.  Memory structures that subserve sentence comprehension.  _Journal of Memory and Language_ , 48(1):67–91.
  * Merkx and Frank (2021) Danny Merkx and Stefan L. Frank. 2021.  Human sentence processing: Recurrence or attention? In _Proceedings of the Workshop on Cognitive Modeling and Computational Linguistics_ , pages 12–22.
  * Nesterov (1983) Yurii E. Nesterov. 1983.  A method for solving the convex programming problem with convergence rate $O\bigl{(}\frac{1}{k^{2}}\bigr{)}$.  _Dokl. Akad. Nauk SSSR_ , 269(3):543–547.
  * Oh et al. (2021) Byung-Doh Oh, Christian Clark, and William Schuler. 2021.  Surprisal estimators for human reading times need character models.  In _Proceedings of the Joint Conference of the 59th Annual Meeting of the Association for Computational Linguistics and the 11th International Joint Conference on Natural Language Processing_ , pages 3746–3757.
  * Oh et al. (2022) Byung-Doh Oh, Christian Clark, and William Schuler. 2022.  Comparison of structural parsers and neural language models as surprisal estimators.  _Frontiers in Artificial Intelligence_ , 5(777963).
  * Parker et al. (2009) Robert Parker, David Graff, Junbo Kong, Ke Chen, and Kazuaki Maeda. 2009.  English Gigaword LDC2009T13.
  * Rachev (1985) Svetlozar Todorov Rachev. 1985.  The Monge-Kantorovich mass transference problem and its stochastic applications.  _Theory of Probability and its Applications_ , 29(4):647–676.
  * Radford et al. (2019) Alec Radford, Jeff Wu, Rewon Child, David Luan, Dario Amodei, and Ilya Sutskever. 2019.  Language models are unsupervised multitask learners.  _OpenAI Technical Report_.
  * Rubner et al. (2000) Yossi Rubner, Carlo Tomasi, and Leonidas J. Guibas. 2000.  The Earth Mover’s Distance as a metric for image retrieval.  _International Journal of Computer Vision_ , 40:99–121.
  * Ryu and Lewis (2021) Soo Hyun Ryu and Richard L. Lewis. 2021.  Accounting for agreement phenomena in sentence comprehension with transformer language models: Effects of similarity-based interference on surprisal and attention.  In _Proceedings of the Workshop on Cognitive Modeling and Computational Linguistics_ , pages 61–71.
  * Ryu and Lewis (2022) Soo Hyun Ryu and Richard L. Lewis. 2022.  Using Transformer language model to integrate surprisal, entropy, and working memory retrieval accounts of sentence processing.  In _35th Annual Conference on Human Sentence Processing_.
  * Sennrich et al. (2016) Rico Sennrich, Barry Haddow, and Alexandra Birch. 2016.  Neural machine translation of rare words with subword units.  In _Proceedings of the 54th Annual Meeting of the Association for Computational Linguistics_ , pages 1715–1725.
  * Shain and Schuler (2021) Cory Shain and William Schuler. 2021.  Continuous-Time Deconvolutional Regression for Psycholinguistic Modeling.  _Cognition_ , 215.
  * Shain et al. (2018) Cory Shain, Marten van Schijndel, and William Schuler. 2018.  Deep syntactic annotations for broad-coverage psycholinguistic modeling.  In _Workshop on Linguistic and Neuro-Cognitive Resources (LREC 2018)_.
  * Smith and Levy (2013) Nathaniel J. Smith and Roger Levy. 2013.  The effect of word predictability on reading time is logarithmic.  _Cognition_ , 128:302–319.
  * van Schijndel and Linzen (2021) Marten van Schijndel and Tal Linzen. 2021.  Single-stage prediction models do not explain the magnitude of syntactic disambiguation difficulty.  _Cognitive Science_ , 45(6).
  * Vaswani et al. (2017) Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N. Gomez, Łukasz Kaiser, and Illia Polosukhin. 2017.  Attention is all you need.  In _Advances in Neural Information Processing Systems_ , volume 30.
  * Wilcox et al. (2020) Ethan Gotlieb Wilcox, Jon Gauthier, Jennifer Hu, Peng Qian, and Roger P. Levy. 2020\.  On the predictive power of neural language models for human real-time comprehension behavior.  In _Proceedings of the 42nd Annual Meeting of the Cognitive Science Society_ , pages 1707–1713.

##  Appendix A CDR Implementation Details

The continuous-time deconvolutional regression (CDR) models used in this work
were fitted using variational inference to estimate the means and variances of
independent normal posterior distributions over all model parameters assuming
an improper uniform prior. Convolved predictors used the three-parameter
ShiftedGamma impuse response function (IRF) kernel:

$f(x;\alpha,\beta,\delta)=\frac{\beta^{\alpha}(x-\delta)^{\alpha-1}e^{-\beta(x-\delta)}}{\Gamma(\alpha)}$
(14)

Posterior means for the IRF parameters were initialized at $\alpha=0.2$,
$\beta=0.5$, and $\delta=-0.2$, which defines a decreasing IRF with a peak
centered at $t=0$ that decays to near-zero within about 1 s. Models were
fitted using the Adam optimizer (Kingma and Ba, 2015) with Nesterov momentum
(Nesterov, 1983; Dozat, 2016), a constant learning rate of 0.001, and
minibatches of size 1,024. For computational efficiency, histories were
truncated at 128 timesteps. Prediction from the network used an exponential
moving average of parameter iterates with a decay rate of 0.999, and the
models were evaluated using maximum a posteriori estimates obtained by setting
all IRF parameters to their posterior means.

For the baseline regression predictors, the ‘index of word position within
each sentence’ predictors were scaled, and the ‘word length in characters’ and
‘saccade length’ predictors were both centered and scaled.


# Model studies of topological phase transitions in materials
with two types of magnetic atoms

Zhuoran He  Wuhan National High Magnetic Field Center $\&$ School of Physics,
Huazhong University of Science and Technology, Wuhan 430074, China  Gang Xu
<EMAIL_ADDRESS>Wuhan National High Magnetic Field Center $\&$ School of
Physics, Huazhong University of Science and Technology, Wuhan 430074, China

###### Abstract

We study the topological phase transitions induced by Coulomb engineering in
three triangular-lattice Hubbard models $AB_{2}$, $AC_{3}$ and $B_{2}C_{3}$,
each of which consists of two types of magnetic atoms with opposite magnetic
moments. The energy bands are calculated using the Schwinger boson method. We
find that a topological phase transition can be triggered by the second-order
(three-site) virtual processes between the two types of magnetic atoms, the
strengths of which are controlled by the on-site Coulomb interaction $U$. This
new class of topological phase transitions have been rarely studied and may be
realized in a variety of real magnetic materials.

##  I Introduction

Topological phase transitions [1, 2, 3, 4] play a key role in condensed matter
physics. Especially, magnetic topological systems [5, 6, 7] often exhibit rich
topological phases due the complicated interplay between electron-electron
interactions, magnetic moments and spin-orbit coupling, which have been
attracting intensive research interests for years [8, 9, 10]. A general model
for the description of magnetic topological insulators is the spin-orbit
coupled Hubbard model [11, 12] with on-site Coulomb interaction $U$. Previous
works on Coulomb engineering and correlation-driven effects in magnetic
topological systems have studied various aspects of this topic including the
Hartree-Fock mean-field theory [6, 13], dynamical screening effects [14], and
phase transitions due to the magnetic exchange coupling [15, 16, 17] using the
Schwinger boson method. These works mostly focus on systems with one type of
magnetic atom, while the topological phase transitions in systems with two
types magnetic atoms are comparatively less studied.

In this paper, we study systems with two types of magnetic atoms [18, 19, 20,
21] with opposite magnetic moments. In such systems, the two types of magnetic
atoms form two sets of Chern bands separately, which then interact via a type
of second-order virtual process of order $\mathcal{O}(t_{1}t_{2}/U)$. These
processes involve the hopping from one type of magnetic atom $i$ to atom $j$
via the other type of magnetic atom $k$ as an intermediate site. We call these
$1/U$-controlled virtual processes the three-site terms, which can induce
interesting topological phase transitions. We study their effect in a 2D
hexagonal Hubbard model with 3 types of lattice sites $A,B,C$ forming
triangular, honeycomb and Kagome sublattices, respectively. By putting spin-up
and spin-down electrons on two of the three types of lattice sites, we
consider $AB_{2}$, $AC_{3}$ and $B_{2}C_{3}$ models and realize
$1/U$-controlled topological phase transitions as characterized by changes in
the Chern numbers of the spin-up and spin-down bands. Our results demonstrate
the interplay between band topology and correlation effects, and present
Coulomb engineering as a powerful tool to manipulate the topological phases of
matter, with potential applications in various solid-state physical systems.

The rest of the paper is organized as follows. In Sec. II, we give the general
formalism of our downfolding technique in the Schwinger boson representation
and obtain the low-energy effective Hamiltonian containing the three-site
terms. In Sec. III, we apply our formalism to the $AB_{2}$, $AC_{3}$ and
$B_{2}C_{3}$ lattice structures to demonstrate the $1/U$-controlled
topological phase transitions. Section IV is a summary and conclusion, with
discussions of potential materials to realize the topological phase
transitions found in our model studies.

##  II Formalism

Suppose an insulating magnetic material is described by the Hubbard model

$\displaystyle
H=\sum_{ij\alpha\beta}t_{ij}^{\alpha\beta}c_{i\alpha}^{\dagger}c_{j\beta}+U\sum_{i}n_{i\uparrow}n_{i\downarrow},$
(1)

where $i,j$ are the site indices and $\alpha,\beta$ label the spin. In the
large $U$ limit, electrons try to avoid double occupancy and thus each site
becomes spin-polarized so as to form different long-range orders such as
ferromagnetism, antiferromagnetism, ferrimagnetism, etc. [22, 23] In the
Schwinger boson representation [24, 25], the electron operator can be
represented as

$\displaystyle c_{i\sigma}^{\dagger}=b_{i\sigma}^{\dagger}h_{i}+\sigma
d_{i}^{\dagger}b_{i\bar{\sigma}},$ (2)

where $\sigma=\;\uparrow(+1),\downarrow(-1)$ is the spin index, $h_{i}$ and
$d_{i}$ are the fermionic holon and doublon operators,
$b_{i\sigma},b_{i\bar{\sigma}}$ are the Schwinger boson operators, and
$\bar{\sigma}=-\sigma$ is the opposite spin of $\sigma$. By using the
downfolding formula

$\displaystyle
H_{\mathrm{eff}}=PHP-\frac{1}{U}PH\bar{P}HP+\mathcal{O}\\!\left(\frac{1}{U^{2}}\right)\\!,$
(3)

where $P$ is the projection operator into the Hilbert space of no doublons, we
obtain the low-energy effective Hamiltonian of the chargeons

$\displaystyle
H_{\mathrm{eff}}=\sum_{ij}\tilde{t}_{ij}h_{i}h_{j}^{\dagger}=\sum_{ij}\tilde{t}_{ij}f_{i}^{\dagger}f_{j}.$
(4)

A particle-hole transformation has been done from the holons $h_{i}\mapsto
f_{i}^{\dagger}$ to the chargeons, with the effective hopping amplitudes
$\tilde{t}_{ij}$ given by

$\displaystyle\tilde{t}_{ij}=\sum_{\alpha\beta}b_{i\alpha}^{\dagger}\\!\left(t_{ij}^{\alpha\beta}-\frac{1}{U}\sum_{k\gamma\delta}\gamma\delta\,t_{ik}^{\alpha\delta}t_{kj}^{\gamma\beta}\,b_{k\bar{\delta}}^{\dagger}b_{k\bar{\gamma}}\right)\\!b_{j\beta}.$
(5)

The derivation of Eqs. (4)–(5) will be given in Appendix A. For magnetically
ordered systems, the bosonic operators can be viewed as $c$-numbers in the
Bose-Einstein condensation (BEC) approximation [26, 27]. Previous works on
topological phase transitions mostly focus on those transitions induced by
changes of the electronic hopping amplitudes $t_{ij}^{\alpha\beta}$, which may
give rise to gap closing, band inversion [28, 29], etc. Here with Eq. (5), we
can study two more types of topological phase transitions in terms of
$\tilde{t}_{ij}$, i.e., a) those induced by changing the magnetic structure
and b) those induced by $1/U$ that controls the strengths of the three-site
virtual processes. This paper focuses on the latter situation. We consider gap
closing of the chargeon bands induced by the change of Hubbard $U$ without
changing the magnetic structure.

For simplicity, we consider a special case for 2D systems that the bare
hopping $t_{ij}^{\alpha\beta}=t_{ij}^{\alpha}\delta_{\alpha\beta}$ conserves
spin and that the magnetic structure is collinear ferrimagnetic in the $z$
direction. Since the magnetic moments have zero $x,y$ components, and no
double occupancy is allowed in the large $U$ limit, every site can be occupied
by either the spin-up electrons or the spin-down electrons only. In such a
situation, Eq. (5) simplifies to

$\displaystyle\tilde{t}_{ij}=\sum_{\sigma}z_{i\sigma}^{*}z_{j\sigma}\left(t_{ij}^{\sigma}-\frac{1}{U}\sum_{k}t_{ik}^{\sigma}t_{kj}^{\sigma}|z_{k\bar{\sigma}}|^{2}\right)\\!,$
(6)

where the bosonic operators $b_{i\sigma}\mapsto z_{i\sigma}$ have been mapped
to $c$-numbers. Now we have a Hamiltonian with two sets of bands formed by
electrons on spin-up sites and spin-down sites, which interact via the second-
order virtual processes described by the three-site $\mathcal{O}(1/U)$ terms.
In this paper, we use Eq. (6) as our simplified formula. Other magnetically
ordered systems with more complex spin configurations such as noncollinear and
spiral spin structures can be studied using Eq. (5).

##  III Results

To study the topological phase transitions within the framework of Eq. (6), we
construct a 2D lattice structure $AB_{2}C_{3}$ with hexagonal symmetry (see
Fig. 1). The $A$ sites form a triangular lattice with one band, which is
topologically trivial. The $B$ sites form a honeycomb lattice with two bands,
which realize the Haldane model. The $C$ sites form a Kagome lattice with
three bands. We will put opposite magnetic moments on two of the three types
of lattice sites and consider electronic phases in the $AB_{2}$, $AC_{3}$ and
$B_{2}C_{3}$ structures, respectively.

Figure 1: The 2D hexagonal lattice structure $AB_{2}C_{3}$. The $A$ sites form
a triangular lattice, the $B$ sites form a honeycomb lattice, and the $C$
sites form a Kagome lattice, all sharing the same lattice vectors
$\vec{a}_{1}$ and $\vec{a}_{2}$.

###  III.1 The AB2 structure

Figure 2: The Chern bands of the $B$ sites (honeycomb) in the $AB_{2}$ model.
Hopping amplitudes $t_{1}=-0.15$ eV, $t_{2}=(0.06+0.04i)$ eV, $t_{3}=-0.01$
eV, $t=0.8$ eV. Hubbard $U=10$ eV in (a) and $U=4$ eV in (b). The Chern
numbers $C_{1,2}$ indicate a topological phase transition (critical $U=5.3$
eV).

We consider an electronic phase with $N_{\uparrow}=N_{\downarrow}=1$ per unit
cell. In case that the on-site orbital energy of an empty $A$ site is lower
than that of an empty $B$ site, one of the spin species (e.g. the $\downarrow$
electrons) would first singly occupy the $A$ sites. Then the other spin
species (the $\uparrow$ electrons) would not occupy the $A$ sites because of
the Hubbard $U$, but instead occupy the $B$ sites at occupancy $0.5$. When the
SOC is considered, the $B$ sites become gapped and the $\uparrow$ electrons
realize the Haldane model with real nearest-neighbor hopping $t_{1}$ and
complex next-nearest-neighbor hopping $t_{2}$. We also consider a real para-
position hopping $t_{3}$ among the $B$ sites and denote the real nearest
$A$-$B$ site hopping as $t$. Following Eq. (6), the effective hoppings
$\tilde{t}_{1-3}$ are given by

$\displaystyle\tilde{t}_{1}=\frac{1}{2}\left(t_{1}-\frac{2t^{2}}{U}\right)\\!,\quad\tilde{t}_{2,3}=\frac{1}{2}\left(t_{2,3}-\frac{t^{2}}{U}\right)\\!.$
(7)

Here we assume the boson fields $z_{A\downarrow}=1$, $z_{A\uparrow}=0$ on the
$A$ sites and $z_{B\uparrow}=1/\sqrt{2}$, $z_{B\downarrow}=0$ on the $B$
sites. In Eq. (6), when the $i,j$ labels are on the $B$ sites, we have
$\sigma=\;\uparrow$ and thus the $k$ label must be on the $A$ sites, which are
occupied by $\bar{\sigma}=\;\downarrow$, so as to mediate a three-site virtual
process $j\rightarrow k\rightarrow i$. All three hoppings $\tilde{t}_{1-3}$
are renormalized by such three-site virtual processes. Due to the three-site-
enhanced hopping $\tilde{t}_{3}$, the $AB_{2}$ model can now realize beyond-
Haldane phases with occupied-band Chern numbers $\pm 2$.

In terms of the effective hoppings $\tilde{t}_{1-3}$, the spin-up Hamiltonian
(i.e., chargeon Hamiltonian restricted to $B$ sites) in atomic gauge takes the
form

$\displaystyle
H_{B}(\vec{k})=\begin{bmatrix}2\,_{\\!}\mathrm{Re}\,_{\\!}[\tilde{t}_{2}\zeta_{2}^{*}(\vec{k})]&\tilde{t}_{1}\zeta_{1}^{*}(\vec{k})+\tilde{t}_{3}\zeta_{1}(2\vec{k})\vspace{0.5ex}\\\
\tilde{t}_{1}\zeta_{1}(\vec{k})+\tilde{t}_{3}\zeta_{1}^{*}(2\vec{k})&2\,_{\\!}\mathrm{Re}\,_{\\!}[\tilde{t}_{2}\zeta_{2}(\vec{k})]\end{bmatrix}\\!,$
(8)

where the functions $\zeta_{1,2}(\vec{k})$ are given by

$\displaystyle\zeta_{1}(\vec{k})$
$\displaystyle=e^{i\vec{k}\cdot\frac{\vec{a}_{1}-\vec{a}_{2}}{3}}+e^{i\vec{k}\cdot\frac{\vec{a}_{1}+2\vec{a}_{2}}{3}}+e^{-i\vec{k}\cdot\frac{2\vec{a}_{1}+\vec{a}_{2}}{3}},$
(9a) $\displaystyle\zeta_{2}(\vec{k})$
$\displaystyle=e^{i\vec{k}\cdot\vec{a}_{1}}+e^{i\vec{k}\cdot\vec{a}_{2}}+e^{-i\vec{k}\cdot(\vec{a}_{1}+\vec{a}_{2})}.$
(9b)

A topological phase transition can be realized as shown in Fig. 2. In Fig. 2a,
the Hubbard $U=10$ eV is large. The two spin species are clearly separated by
the Hubbard interaction with almost forbidden three-site virtual hoppings. The
spin-up electrons form a Haldane phase on the $B$ sites with occupied-band
Chern number $C_{1}=+1$ and unoccupied-band Chern number $C_{2}=-1$. The spin-
down electrons fully occupy the triangular sites ($A$ sites) and form a
topologically trivial band (not plotted) with Chern number $0$. As $U$ gets
smaller, the three-site virtual processes $\sim\mathcal{O}(1/U)$ become
stronger and the para-position hopping $\tilde{t}_{3}$ is significantly
enhanced. The band gap in Fig. 2a then closes at the M point at critical
$U=5.3$ eV and reopens as $U$ is further reduced to form a beyond-Haldane
phase with $C_{1}=-2$ and $C_{2}=+2$ (see Fig. 2b for $U=4$ eV). Since the
contribution $t^{2}/U$ of the second-order virtual processes is real, the
imaginary part $\mathrm{Im}\,\tilde{t}_{2}=\mathrm{Im}\,t_{2}$ remains
unaffected by $U$. Therefore, the system can undergo topological phase
transitions between $C_{1}=+1\leftrightarrow-2$ (if $\mathrm{Im}\,t_{2}>0$) or
$C_{1}=-1\leftrightarrow+2$ (if $\mathrm{Im}\,t_{2}<0$), but not in between
the $C_{1}=\pm 1$ (or $\pm 2$) phases by tuning $U$.

###  III.2 The AC3 structure

Figure 3: The Chern bands of the $C$ sites (Kagome) in the $AC_{3}$ structure.
Hopping amplitudes $t_{1}=-(0.6+0.2i)$ eV, $t_{2}=(0.1+0.1i)$ eV,
$t_{3}=-0.25$ eV, $t=0.8$ eV. Hubbard $U=10$ eV in (a) and $U=4$ eV in (b).
The Chern numbers $C_{1-3}$ indicate a topological phase transition (critical
$U=4.9$ eV).

Consider an electronic phase in which the $A$ sites are singly occupied by the
$\downarrow$ electrons and the $C$ sites are occupied by the $\uparrow$
electrons at occupancy $1/3$. The situation is similar to $AB_{2}$, except
that the $C$ sites form a Kagome lattice. We consider the nearest-neighbor and
next-nearest-neighbor hoppings $t_{1}$, $t_{2}$, and real para-position
hopping $t_{3}$ of the $C$-site hexagons. Both $t_{1}$ and $t_{2}$ can be
complex. The real nearest-neighbor $A$-$C$ site hopping is denoted as $t$.
From Eq. (6), we have

$\displaystyle\tilde{t}_{1-3}=\frac{1}{3}\left(t_{1-3}-\frac{t^{2}}{U}\right)\\!,$
(10)

assuming $z_{A\downarrow}=1$, $z_{A\uparrow}=0$ for the $A$ sites and
$z_{C\uparrow}=1/\sqrt{3}$, $z_{C\downarrow}=0$ for the $C$ sites. In terms of
$\tilde{t}_{1-3}$, the $C$-site Kagome Hamiltonian takes the form

$\displaystyle H_{C}(\vec{k})=\sum_{\nu=1}^{3}H_{C}^{(\nu)}(\vec{k}),$ (11)

where the nearest-neighbor hopping $H_{C}^{(1)}(\vec{k})$, the next-nearest-
neighbor hopping $H_{C}^{(2)}(\vec{k})$ and the para-position hopping
$H_{C}^{(3)}(\vec{k})$ Hamiltonians are given specifically in Appendix B. A
topological phase transition analogous to the $AB_{2}$ situation is realized
in Fig. 3. In Fig. 3a, the Hubbard $U=10$ eV is large and the three-site
virtual hoppings are almost forbidden. As $U$ gets smaller, the para-position
hopping $\tilde{t}_{3}$ is significantly enhanced. The occupied-band Chern
number changes from $C_{1}=+1$ (see Fig. 3a) to $C_{1}=-2$ (see Fig. 3b) when
the gap closes at the M point at critical $U=4.9$ eV. In the mean time, the
Chern number $C_{2}$ of the middle band changes from $0$ to $+3$ and the Chern
number of the flat band on the top $C_{3}=-1$ remains unchanged.

In Secs. III.1–III.2, we have studied the enhancement effect of para-position
hopping $\tilde{t}_{3}$ due the three-site virtual processes proportional to
$1/U$. We find that in both the honeycomb and the Kagome lattices, the three-
site processes can lead to topological phase transitions of
$C_{1}=+1\leftrightarrow-2$ (or symmetrically $C_{1}=-1\leftrightarrow+2$) by
closing the band gap at the M point. Because $t^{2}/U$ is real, we cannot
realize topological phase transitions between the $C_{1}=\pm 1$ phases. We
will demonstrate in Sec. III.3 that the $+1\leftrightarrow-1$ transitions can
be realized in the B2C3 model by making the contributions of the three-site
processes $\mathcal{O}(tt^{\prime}/U)$ complex.

###  III.3 The B2C3 structure

Figure 4: The hoppings considered in B2C3. Here $t,t^{\prime}$ are between the
$B$-$C$ sites and $t_{1},t_{2}$ are the nearest-neighbor and next-nearest-
neighbor hoppings of the $C$-site Kagome lattice. All hoppings except $t$ can
be complex due to the SOC. The para-position hoppings are ignored.

In this section, we consider an electronic phase with
$N_{\uparrow}=N_{\downarrow}=2$ per unit cell. Let the two $B$ sites in a unit
cell be singly occupied by the $\downarrow$ electrons and the three $C$ sites
be occupied by $\uparrow$ electrons at occupancy $2/3$. We consider the
hoppings $t,t^{\prime}$ between the $B$-$C$ sites and hoppings $t_{1},t_{2}$
among the $C$ sites as shown in Fig. 4. Since the total Chern number of the
two spin-down bands on the $B$ sites is zero (c.f. Fig. 2a), we focus on the
topological properties of the Kagome bands, which are controlled by $1/U$.
From Eq. (6), we have

$\displaystyle\tilde{t}_{1}=\frac{2}{3}\left(t_{1}-\frac{t^{2}}{U}\right)\\!,\quad\tilde{t}_{2}=\frac{2}{3}\left(t_{2}-\frac{2tt^{\prime}}{U}\right)\\!,$
(12)

assuming $z_{B\downarrow}=1$, $z_{B\uparrow}=0$ and
$z_{C\uparrow}=\sqrt{2/3}$, $z_{C\downarrow}=0$. The para-position hopping
$\tilde{t}_{3}=0$ of the Kagome lattice is ignored. Even though
$\tilde{t}_{3}$ can be mediated by $t^{\prime 2}/U$, these contributions are
small assuming $|t|\gg|t^{\prime}|$. Notice that the $B$-$C$ site hopping $t$
is real, while $t^{\prime}$ can be complex due to SOC. We define for
$\uparrow$ electrons that the blue-line hoppings in Fig. 4 is $t^{\prime}$ in
clockwise directions and $(t^{\prime})^{*}$ in counter-clockwise directions.
The Hamiltonian $H_{C}(\vec{k})$ is still given in Appendix B with the
effective hoppings $\tilde{t}_{1,2}$ now given by Eq. (12). A topological
phase transition is realized as shown in Fig. 5.

Figure 5: The Chern bands of the $C$ sites in the $B_{2}C_{3}$ model. Hopping
amplitudes $t_{1}=(0.6-0.1i)$ eV, $t_{2}=-(0.1+0.02i)$ eV, $t=0.8$ eV,
$t^{\prime}=(0.1+0.1i)$ eV. Hubbard $U=10$ eV in (a) and $U=4$ eV in (b). The
Chern numbers $C_{1-3}$ indicate a topological phase transition (critical
$U=5.3$ eV).

In Fig. 5a, the Hubbard $U=10$ eV and the occupied-band Chern number
$C_{1}+C_{2}=+1$, which is determined by the imaginary parts
$\mathrm{Im}\,\tilde{t}_{1,2}$ of the effective hoppings in the Kagome
lattice. As $U$ gets smaller, since $t^{2}/U$ is real,
$\mathrm{Im}\,\tilde{t}_{1}$ remains unchanged, so only the $tt^{\prime}/U$
term in Eq. (12) can affect $\mathrm{Im}\,\tilde{t}_{2}$. The band gap closes
at the K point at critical $U=5.3$ eV and then reopens to give rise to a
$C_{1}+C_{2}=-1$ phase as $U$ further decreases to $4$ eV (see Fig. 5b). The
Chern number of the flat band at the bottom $C_{1}=+1$ remains unchanged
throughout the process. Because the imaginary part of the hopping amplitudes
can be tuned by $1/U$, the phase separation between Chern numbers $\pm 1$ is
broken. A topological phase transition between the $\pm 1$ phases can now be
realized by tuning the Hubbard $U$ due to the complex virtual hopping
$\mathcal{O}(tt^{\prime}/U)$.

##  IV Conclusion

We have demonstrated in this paper that the three-site virtual processes in
the large $U$ limit of the Hubbard model can exhibit interesting
renormalization effects of the hopping amplitudes and give rise to topological
phase transitions in the low-energy effective theory. We constructed 2D
lattice models to realize the $1/U$ control of the honeycomb and Kagome
lattices. In the $AB_{2}$ model, a topological phase transition between the
Haldane phase and beyond-Haldane phase is realized by considering the
enhancement effect of the para-position hopping $\tilde{t}_{3}$ due to the
$A$-site mediated virtual hoppings proportional to $1/U$. The $AC_{3}$ model
realizes a similar phase transition on the Kagome lattice. Both transitions
close the band gap at the M point. In the $B_{2}C_{3}$ model, we also realize
topological phase transitions on the Kagome lattice, but the band gap closes
at the K point. The contribution $\mathcal{O}(tt^{\prime}/U)$ of the three-
site processes can be complex and drives the system across the phases boundary
of occupied-band Chern number $=\pm 1$.

The phase transitions found in our model studies are realized using collinear
antiferromagnetic (or ferrimagnetic) spin configurations. The spin-up and
spin-down electrons occupy inequivalent lattice sites. In the examples shown
in this paper, for simplicity, we let one spin species fully occupy one type
of lattice site so as to be topologically trivial, and use them to control the
topological phase of the other spin species via the three-site terms.
Interesting directions for further studies could be that both spin species
exhibit topological properties and mutually influence each other via the
three-site terms, or the realization of similar Coulomb engineering effects in
non-collinear spin systems.

Figure 6: Possible realizations of the $AB_{2}C_{3}$ lattice in a 3D hexagonal
crystal structure with (a) alternating $AB_{2}$ and $C_{3}$ layers and (b)
alternating $AC_{3}$ and $B_{2}$ layers. Both structures have the P6/mmm space
group symmetry.

Finally, we would like to discuss the possible realizations of our model in
real materials. The $1/U$-controlled topological phase transitions can be
realized without restricting the atoms to the same 2D plane. Two possible 3D
structures are shown in Fig. 6, both with P6/mmm symmetry. Examples of
materials with the structure of Fig. 6a are $R$Co3B2 [30] with $R=\,$rare-
earth elements, and also GdNi3Ga2 [31], etc, which are potential candidates
for the $AC_{3}$ model of two types of magnetic atoms. Coplanar $AC_{3}$
candidates include the TiNi3-type compounds [32, 33, 34] with shifted layers
of close-packed $AC_{3}$ structures. Candidates for the $AB_{2}$ model include
e.g. UNi2Al3 [35] and EuCo2Al9 [36], etc. We expect our work to be interesting
to the fields of magnetism in alloys, ferrimagnets and other materials with
multiple types of magnetic atoms.

## Acknowledgments

This work is supported by the National Key Research and Development Program of
China (2018YFA0307000), and the National Natural Science Foundation of China
(11874022). Z. H would thank the support of the 66th Chinese Postdoc
Fellowship. We would also like to thank the helpful discussions with Prof.
Biao Lian at Princeton Center for Theoretical Science in Princeton University
at the early stage of this work.

##  Appendix A Derivation of the low-energy effective Hamiltonian

By plugging Eq. (1) into Eq. (3), one obtains

$\displaystyle H_{\mathrm{eff}}$
$\displaystyle=\sum_{ij\alpha\beta}t_{ij}^{\alpha\beta}Pc_{i\alpha}^{\dagger}c_{j\beta}P-\frac{1}{U}\sum_{ijkl}\sum_{\alpha\beta\gamma\delta}t_{ij}^{\alpha\beta}t_{kl}^{\gamma\delta}Pc_{i\alpha}^{\dagger}c_{j\beta}\bar{P}c_{k\gamma}^{\dagger}c_{l\delta}P.$
(13)

Then plugging in Eq. (2), one finds that the projections $P$ pick out the
following terms

$\displaystyle
H_{\mathrm{eff}}=\sum_{ij\alpha\beta}t_{ij}^{\alpha\beta}h_{i}h_{j}^{\dagger}b_{i\alpha}^{\dagger}b_{j\beta}-\frac{1}{U}\sum_{ijkl}\sum_{\alpha\beta\gamma\delta}\beta\gamma\,t_{ij}^{\alpha\beta}t_{kl}^{\gamma\delta}\,b_{i\alpha}^{\dagger}b_{j\bar{\beta}}^{\dagger}\,b_{k\bar{\gamma}}\,b_{l\delta}\,h_{i}\,(Pd_{j}d_{k}^{\dagger}P)\,h_{l}^{\dagger}.$
(14)

The bosonic operators are automatically normal ordered. Notice that
$Pd_{j}d_{k}^{\dagger}P=\delta_{jk}$, because the doublon created must also be
the doublon destructed so as to go back to the no-doublon subspace. One may
then set $j=k$ and rename the dummy indices $l\mapsto j$ and
$\beta\leftrightarrow\delta$ to obtain

$\displaystyle
H_{\mathrm{eff}}=\sum_{ij}h_{i}h_{j}^{\dagger}\left[\sum_{\alpha\beta}b_{i\alpha}^{\dagger}\\!\left(t_{ij}^{\alpha\beta}-\frac{1}{U}\sum_{k}\sum_{\gamma\delta}\gamma\delta\,t_{ik}^{\alpha\delta}t_{kj}^{\gamma\beta}\,b_{k\bar{\delta}}^{\dagger}\,b_{k\bar{\gamma}}\right)\\!b_{j\beta}\right]\\!.$
(15)

This result agrees with Eqs. (4)–(5) in the main text by defining the quantity
in the square bracket as $\tilde{t}_{ij}$. All $\mathcal{O}(1/U)$
renormalizations of $\tilde{t}_{ij}$ are considered in this formalism. Then we
do a particle-hole transformation $h_{i}\mapsto f_{i}^{\dagger}$ to the holon
operators and map the bosonic operators $b_{i\sigma}\mapsto z_{i\sigma}$ to
$c$-numbers and obtain

$\displaystyle
H_{\mathrm{eff}}=\sum_{ij}f_{i}^{\dagger}f_{j}\left[\sum_{\alpha\beta}z_{i\alpha}^{*}z_{j\beta}\left(t_{ij}^{\alpha\beta}-\frac{1}{U}\sum_{k\gamma\delta}\gamma\delta\,t_{ik}^{\alpha\delta}t_{kj}^{\gamma\beta}\,{}_{\\!}z_{k\bar{\delta}}^{*}z_{k\bar{\gamma}}\right)\right]=\sum_{ij}\tilde{t}_{ij}f_{i}^{\dagger}f_{j}.$
(16)

In the special case that the bare hopping
$t_{ij}^{\alpha\beta}=t_{ij}^{\alpha}\delta_{\alpha\beta}$ conserves spin, we
have

$\displaystyle\tilde{t}_{ij}=\sum_{\alpha\beta}z_{i\alpha}^{*}z_{j\beta}\left(t_{ij}^{\alpha}\delta_{\alpha\beta}-\frac{1}{U}\sum_{k}\alpha\beta\,t_{ik}^{\alpha}t_{kj}^{\beta}z_{k\bar{\alpha}}^{*}z_{k\bar{\beta}}\right)\\!.$
(17)

Then the collinear ferrimagnetic structure in the $z$ direction (perpendicular
to the 2D lattice plane) with no double occupancy eliminates the
$\alpha\neq\beta$ terms because site $k$ can only be occupied by one type of
spin species. Therefore, one obtains Eq. (6) in the main text.

##  Appendix B Kagome Hamiltonian in terms of $\tilde{t}_{1-3}$

In terms of the effective hoppings $\tilde{t}_{1-3}$, the full Kagome
Hamiltonian $H_{C}(\vec{k})$ contains 3 parts as defined by Eq. (11): the
nearest-neighbor hopping Hamiltonian as given by

$\displaystyle
H_{C}^{(1)}(\vec{k})=\begin{bmatrix}0&2\tilde{t}_{1}\cos\left(\vec{k}\cdot\frac{\vec{a}_{1}+\vec{a}_{2}}{2}\right)&2\tilde{t}_{1}^{*}\cos\left(\vec{k}\cdot\frac{\vec{a}_{2}}{2}\right)\\\
2\tilde{t}_{1}^{*}\cos\left(\vec{k}\cdot\frac{\vec{a}_{1}+\vec{a}_{2}}{2}\right)&0&2\tilde{t}_{1}\cos\left(\vec{k}\cdot\frac{\vec{a}_{1}}{2}\right)\\\
2\tilde{t}_{1}\cos\left(\vec{k}\cdot\frac{\vec{a}_{2}}{2}\right)&2\tilde{t}_{1}^{*}\cos\left(\vec{k}\cdot\frac{\vec{a}_{1}}{2}\right)&0\end{bmatrix}\\!,$
(18)

the next-nearest-neighbor hopping Hamiltonian as given by

$\displaystyle
H_{C}^{(2)}(\vec{k})=\begin{bmatrix}0&2\tilde{t}_{2}\cos\left(\vec{k}\cdot\frac{\vec{a}_{1}-\vec{a}_{2}}{2}\right)&2\tilde{t}_{2}^{*}\cos\left[\vec{k}\cdot\left(\vec{a}_{1}+\frac{\vec{a}_{2}}{2}\right)\right]\\\
2\tilde{t}_{2}^{*}\cos\left(\vec{k}\cdot\frac{\vec{a}_{1}-\vec{a}_{2}}{2}\right)&0&2\tilde{t}_{2}\cos\left[\vec{k}\cdot\left(\frac{\vec{a}_{1}}{2}+\vec{a}_{2}\right)\right]\\\
2\tilde{t}_{2}\cos\left[\vec{k}\cdot\left(\vec{a}_{1}+\frac{\vec{a}_{2}}{2}\right)\right]&2\tilde{t}_{2}^{*}\cos\left[\vec{k}\cdot\left(\frac{\vec{a}_{1}}{2}+\vec{a}_{2}\right)\right]&0\end{bmatrix}\\!,$
(19)

and the para-position hopping Hamiltonian as given by

$\displaystyle
H_{C}^{(3)}(\vec{k})=\begin{bmatrix}2\tilde{t}_{3}\cos\left(\vec{k}\cdot\vec{a}_{1}\right)&0&0\\\
0&2\tilde{t}_{3}\cos\left(\vec{k}\cdot\vec{a}_{2}\right)&0\\\
0&0&2\tilde{t}_{3}\cos\left[\vec{k}\cdot(\vec{a}_{1}+\vec{a}_{2})\right]\end{bmatrix}\\!,$
(20)

all written in the atomic gauge. In the $AC_{3}$ model, we consider the
$1/U$-control of all three effective hoppings $\tilde{t}_{1-3}$. In the B2C3
model, we restrict ourselves to $\tilde{t}_{1-2}$, i.e., setting
$\tilde{t}_{3}=0$.

## References

  * Sen  _et al._ [2020] S. Sen, P. J. Wong, and A. K. Mitchell, The mott transition as a topological phase transition, Phys. Rev. B 102, 081110 (2020).
  * Rufo  _et al._ [2019] S. Rufo, N. Lopes, M. A. Continentino, and M. A. R. Griffith, Multicritical behavior in topological phase transitions, Phys. Rev. B 100, 195432 (2019).
  * Mitchell  _et al._ [2018] N. P. Mitchell, L. M. Nash, and W. T. M. Irvine, Realization of a topological phase transition in a gyroscopic lattice, Phys. Rev. B 97, 100302 (2018).
  * Wang and Fu [2017] Y. Wang and L. Fu, Topological phase transitions in multicomponent superconductors, Phys. Rev. Lett. 119, 187003 (2017).
  * Nevola  _et al._ [2020] D. Nevola, H. X. Li, J.-Q. Yan, R. G. Moore, H.-N. Lee, H. Miao, and P. D. Johnson, Coexistence of surface ferromagnetism and a gapless topological state in ${\mathrm{mnbi}}_{2}{\mathrm{te}}_{4}$, Phys. Rev. Lett. 125, 117205 (2020).
  * Ivanov  _et al._ [2019] V. Ivanov, X. Wan, and S. Y. Savrasov, Topological insulator-to-weyl semimetal transition in strongly correlated actinide system unisn, Phys. Rev. X 9, 041055 (2019).
  * Wang  _et al._ [2015a] J. Wang, B. Lian, and S.-C. Zhang, Electrically tunable magnetism in magnetic topological insulators, Phys. Rev. Lett. 115, 036805 (2015a).
  * Zhai and Jin [2014] X. Zhai and G. Jin, Photoinduced topological phase transition in epitaxial graphene, Phys. Rev. B 89, 235416 (2014).
  * Wang  _et al._ [2014] L. Wang, H.-H. Hung, and M. Troyer, Topological phase transition in the hofstadter-hubbard model, Phys. Rev. B 90, 205111 (2014).
  * Cai  _et al._ [2008] Z. Cai, S. Chen, S. Kou, and Y. Wang, Properties of a class of topological phase transitions, Phys. Rev. B 78, 035123 (2008).
  * De Franco  _et al._ [2018] C. De Franco, L. F. Tocchio, and F. Becca, Metal-insulator transitions, superconductivity, and magnetism in the two-band hubbard model, Phys. Rev. B 98, 075117 (2018).
  * Misumi  _et al._ [2017] K. Misumi, T. Kaneko, and Y. Ohta, Mott transition and magnetism of the triangular-lattice hubbard model with next-nearest-neighbor hopping, Phys. Rev. B 95, 075124 (2017).
  * Zhu  _et al._ [2019] Q. Zhu, Q. Tong, H. Sun, Y. Wang, and W. Yao, Coulomb effects on topological band inversion in the moiré of WSe2/BAs heterobilayer, 2D Materials 6, 045037 (2019).
  * Rösner and Lado [2020] M. Rösner and J. L. Lado, Coulomb-engineered topology (2020), arXiv:2008.07990 [cond-mat.str-el] .
  * Lei and Jun [2007] X. Lei and Z. Jun, Slave-fermion mean-field theory of heisenberg model, Communications in Theoretical Physics 48, 363 (2007).
  * Yoshioka [1989] D. Yoshioka, Slave-fermion mean field theory of the hubbard model, Journal of the Physical Society of Japan 58, 1516 (1989).
  * Jayaprakash  _et al._ [1989] C. Jayaprakash, H. R. Krishnamurthy, and S. Sarker, Mean-field theory for the t-j model, Phys. Rev. B 40, 2610 (1989).
  * Ohnishi  _et al._ [1983] H. Ohnishi, K. Katoh, and K. Motizuki, Magnetic phase diagram of system with two different atoms in unit cell, Journal of Magnetism and Magnetic Materials 31-34, 55 (1983).
  * MARCOS  _et al._ [1990] M. MARCOS, P. ROMERO, and J.-L. SERRANO, Cheminform abstract: Nematic liquid crystal materials containing nickel(ii) or copper(ii) atoms. two different kinds of magnetic behaviour., ChemInform 21 (1990), https://onlinelibrary.wiley.com/doi/pdf/10.1002/chin.199016323 .
  * Wang  _et al._ [2015b] P. Wang, W. Ren, L. Bellaiche, and H. Xiang, Predicting a ferrimagnetic phase of zn2feoso6 with strong magnetoelectric coupling, Phys. Rev. Lett. 114, 147204 (2015b).
  * Halder  _et al._ [2010] M. Halder, S. M. Yusuf, M. D. Mukadam, and K. Shashikala, Magnetocaloric effect and critical behavior near the paramagnetic to ferrimagnetic phase transition temperature in ${\text{tbco}}_{2-x}{\text{fe}}_{x}$, Phys. Rev. B 81, 174402 (2010).
  * Gouveia and Dias [2015] J. Gouveia and R. Dias, Magnetic phase diagram of the hubbard model in the lieb lattice, Journal of Magnetism and Magnetic Materials 382, 312 (2015).
  * Wel  _et al._ [1989] G. Z. Wel, H. Q. Nie, and K. Y. Zhang, Correlation effects of the antiferromagnetic phase in the hubbard model, physica status solidi (b) 156, 363 (1989).
  * Wang and Tao [2000] F. Wang and R. Tao, Schwinger-boson approach to the two-dimensional antiferromagnetic heisenberg model beyond the mean-field approximation, Phys. Rev. B 61, 3508 (2000).
  * Rech  _et al._ [2006] J. Rech, P. Coleman, G. Zarand, and O. Parcollet, Schwinger boson approach to the fully screened kondo model, Phys. Rev. Lett. 96, 016601 (2006).
  * Ueda and Totsuka [2007] H. T. Ueda and K. Totsuka, Ground-state phase diagram and magnetic properties of a tetramerized spin-$1/2$ ${J}_{1}-{J}_{2}$ model: Bec of bound magnons and absence of the transverse magnetization, Phys. Rev. B 76, 214428 (2007).
  * Davis  _et al._ [2016] M. J. Davis, T. M. Wright, T. Gasenzer, S. A. Gardiner, and N. P. Proukakis, Formation of bose-einstein condensates (2016), arXiv:1601.06197 [cond-mat.quant-gas] .
  * Wojek  _et al._ [2014] B. M. Wojek, P. Dziawa, B. J. Kowalski, A. Szczerbakow, A. M. Black-Schaffer, M. H. Berntsen, T. Balasubramanian, T. Story, and O. Tjernberg, Band inversion and the topological phase transition in (pb,sn)se, Phys. Rev. B 90, 161202 (2014).
  * Wang  _et al._ [2017] D. Wang, L. Chen, H. Liu, C. Shi, X. Wang, G. Cui, P. Zhang, and Y. Chen, Strain induced band inversion and topological phase transition in methyl-decorated stanene film, Scientific Reports 7, 17089 (2017).
  * Ido  _et al._ [1999] H. Ido, O. Nashima, T. Ito, T. Kaneda, Y. Saito, K. Konno, H. Yoshida, and M. Motokawa, Magnetic properties of rco3b2 (r=rare-earth) compounds, Journal of Applied Physics 85, 4865 (1999), https://doi.org/10.1063/1.370047 .
  * Pechev  _et al._ [2000] S. Pechev, J.-L. Bobet, B. Chevalier, B. Darriet, and F. Weill, Relationship between the synthesis conditions and the structural, electrical, magnetic, and hydrogen absorption properties of the ternary compounds gdni3x2 with x=al, ga, or sn, Journal of Solid State Chemistry 150, 62 (2000).
  * Xu  _et al._ [1993] J.-h. Xu, W. Lin, and A. J. Freeman, Electronic structure and phase stability of ${\mathit{a}}_{3}$ti (a=fe, co, ni, and cu), Phys. Rev. B 48, 4276 (1993).
  * Wodniecka  _et al._ [1995] B. Wodniecka, P. Wodniecki, M. Marszalek, and A. Hrynkiewicz, Electric quadrupole interaction at 181ta probes in hexagonal tini3-type compounds, Journal of Alloys and Compounds 219, 128 (1995), eleventh international conference on solid compounds of transition elements.
  * Le  _et al._ [2010] M. D. Le, K. A. McEwen, E. Colineau, J.-C. Griveau, and R. Eloirdi, Magnetic and electrical properties of $(\text{Pu},\text{Lu}){\text{pd}}_{3}$, Phys. Rev. B 82, 155136 (2010).
  * Geibel  _et al._ [1993] C. Geibel, A. Böhn, R. Caspary, K. Gloos, A. Grauel, P. Hellmann, R. Modler, C. Schank, G. Weber, and F. Steglich, Ground state properties of uni2al3 and upd2al3, Physica B: Condensed Matter 186-188, 188 (1993).
  * Thiede and Jeitschko [1999] V. . T. Thiede and W. Jeitschko, Crystal structure of europium cobalt aluminide (1/2/9), euco2al9, Zeitschrift für Kristallographie - New Crystal Structures 214, 149 (1999).


# MIDI Miner – A Python Library for tonal tension and track classification

We present a Python library, called Midi Miner, that can calculate tonal
tension and classify different tracks. MIDI (Music Instrument Digital
Interface) is a hardware and software standard for communicating musical
events between digital music devices. It is often used for tasks such as music
representation, communication between devices, and even music generation [5].
Tension is an essential element of the music listening experience, which can
come from a number of musical features including timbre, loudness and harmony
[3]. Midi Miner provides a Python implementation for the tonal tension model
based on the spiral array [1] as presented by Herremans and Chew [4]. Midi
Miner also performs key estimation and includes a track classifier that can
disentangle melody, bass, and harmony tracks. Even though tracks are often
separated in MIDI files, the musical function of each track is not always
clear. The track classifier keeps the identified tracks and discards messy
tracks, which can enable further analysis and training tasks.

The spiral array [1] is a three-dimensional mathematical model of tonality
developed by Prof. Elaine Chew. It consists of a helix for pitch classes,
chords, and keys. It is important to note that this model takes pitch spelling
into account, meaning that Ab and G# will have a different position in the
pitch class helix. Since MIDI files do not capture pitch spelling, we include
an algorithm for estimating pitch spelling. First, the key index of the MIDI
file is set as the most frequently occurring pitch class. The mode of the key
(major or minor) cannot be differentiated at this stage, but the key index can
inform the pitch spelling. Based on the estimated key index, the piece will
fall into one of two classes, the first is based on the key indices C, D, E,
G, A, and B, and the second on Db, Eb, F, Gb, Ab ad Bb. The accidentals of
pieces in the first class will be mapped to sharps if needed, whereas the
accidentals of pieces in the second class will be mapped to flats. After the
pitch spelling has been estimated for all note, the final key can be
determined using the spiral array key detection algorithm [2].

Tonal tension [4] define three measures of tonal tension: cloud diameter
(dissonance), cloud momentu (tonal movement), and tensile strain (distance to
the key). Mini Miner implements these three tension measures based on a
default window length of one half note (customizable by the user). The
implemented model works for polyphonic music, by weighing the occurrence of
each pitch class per time window. After inputting a MIDI file, Midi Miner also
provides a graphical representation of the three tension measures. Figure 1
provides an example output for the tensile strain measure (i.e., distance to
the key). A sliding window of 16 beats with average tensile strain is used to
detect potential key changes by comparing the ratio between the mean tensile
strain and the adjacent window. If this $ratio>2$ for four consecutive
windows, then a key change is flagged. In Figure 1, Midi Miner detects a key
change in bar 61.

Figure 1: Tensile strain for each bar of the MIDI file available
online111https://github.com/ruiguo-bio/midi-miner/blob/master/example/new.mid.

Midi track classifier Midi Miner also implements a track classifier that
untangles melody, bass, and harmony tracks. These tracks are of interest for
anybody preparing a dataset for, e.g., music generation. A random forest
classifier was trained by using 5,006 Chinese pop songs collected from the
Internet. Among those files, 3,265 files have at least one track label. These
MIDI track labels are strings marked by the author, which usually provide
information about the track type. By looking at the labels in those MIDI
files, we manually added specific labels for melody, bass, chords, which were
used to train our classifier. The drum tracks are currently discarded. Similar
to [6], who select 34 features to classify the melody track, we computed 30
features for each track. These are used as input to train our classifier and
are described in detail on our github page. The performance of the track
classifier is evaluated in Table 1 using a test set (25% or the original
dataset). Based on the classification result, Midi Miner will extract the
melody, bass and harmony tracks and output a new MIDI file with only those
tracks.

track | prediction | precision | recall | F1 score
---|---|---|---|---
melody | False | 0.99 | 0.99 | 0.99
| True | 0.95 | 0.93 | 0.94
bass | False | 0.98 | 1.00 | 0.99
| True | 0.97 | 0.92 | 0.95
harmony | False | 0.98 | 0.98 | 0.98
| True | 0.92 | 0.89 | 0.91
Table 1: Track classification result

Download Midi Miner The Midi Miner library is available
online222https://github.com/ruiguo-bio/midi-miner. In the future, we aim to
include a more comprehensive extraction of MIDI features, as well as more
complete music analysis tools.

ACKNOWLEDGMENTS

This work has received funding from the Chinese scholarship council and MOE T2
grant number MOE2018-T2-2-161.

## References

  * [1] E. Chew.  Mathematical and computational modeling of tonality.  AMC, 10:12, 2014.
  * [2] C.-H. Chuan and E. Chew.  Polyphonic audio key finding using the spiral array ceg algorithm.  In 2005 IEEE International Conference on Multimedia and Expo, pages 21–24. IEEE, 2005.
  * [3] M.M. Farbood.  A parametric, temporal model of musical tension.  Music Perception, 29(4):387–428, 2012.
  * [4] D. Herremans and E. Chew.  Tension ribbons: Quantifying and visualising tonal tension.  In Second International Conference on Technologies for Music Notation and Representation (TENOR), volume 2, pages 8–18, Cambridge, UK, 05/2016 2016.
  * [5] D. Herremans, C.-H. Chuan, and E. Chew.  A functional taxonomy of music generation systems.  ACM Computing Surveys (CSUR), 50(5):69, 2017.
  * [6] David Rizo, Pedro J Ponce De León, Antonio Pertusa, Carlos Pérez-Sancho, and José Manuel Iñesta Quereda.  Melody track identification in music symbolic files.  In FLAIRS Conference, pages 254–259, 2006.