#  Study of quantum non-locality by CHSH function and its extension
in disordered fermions

Yoshihito Kuno1 1Graduate School of Engineering Science, Akita University,
Akita 010-8502, Japan

###### Abstract

Quantum non-locality is an important concept in quantum physics. In this work,
we study the quantum non-locality in a fermion many-body system under quasi-
periodic disorders. The Clauser-Horne-Shimony-Holt (CHSH) inequality is
systematically investigated, which quantifies quantum non-locality between two
sites. We find that the quantum non-locality explicitly characterize the
extended and critical phase transitions, and further that in the globally
averaged picture of maximum value of the quantum non-locality the CHSH
inequality is not broken, but for a local pair in the internal of the system
the violation probability of the CHSH inequality becomes sufficiently finite.
Further we investigate an extension of the CHSH inequality, Mermin-Klyshko-
Svetlichny (MKS) polynomials, which can characterize multipartite quantum non-
locality. We also find a similar behavior to the case of CHSH inequality. In
particular, in the critical regime and on a transition point, the adjacent
three qubit MKS polynomial in a portion of the system exhibits a quantum non-
local violation regime with a finite probability.

##  I Introduction

Study of quantum correlation is one of the most attracted research field in
physics Horodecki_family . In particular, the characterization of quantum
correlation for many body system is now on-going problem. Quantum entanglement
is one of the most famous measure of quantum correlation in many-body system.
The behaviors of quantum entanglement or entanglement entropy have been
extensively studied Horodecki_family ; Eisert2010 . They succeeded in
characterizing non-trivial phases such as symmetry protected topological
phase, topological ordered phase Wen_text .

There is another concept elucidating quantum correlation. It is quantum non-
locality. Bell inequality Bell1969  is the first celebrated example of it by
exhibiting a violation of its inequality. Further, the simpler version of it
called Clauser, Horne, Shimony, and Holt (CHSH) inequality CHSH1969 ,
elucidating quantum non-locality between two spins the measurement outcome of
which is dichotomic. For a singlet entangled state, the CHSH inequality
explicitly breaks, describing a boundary line between quantum non-locality and
local classical phenomena described by a hidden variable theory. Here we are
confronted with two concepts of correlation, quantum entanglement and quantum
non-locality. An important issue is to clarify how these physical concepts are
related or different from each other. Entanglement may not be identical to the
concept of quantum non-locality measured by Bell or CHSH inequality. As an
amazing example, a mixed state with entanglement does not exhibit quantum non-
locality Werner1989 ; Augusiak2015 , while entanglement and quantum non-
locality can be identical in two qubit system Gisin1991 . Further study about
this issue is required.

The study of quantum non-locality in many-body systems has been less-attracted
than that of many-body entanglement. In particular, in a highly quantum
correlated state, how quantum non-locality and further multipartite one appear
and what structure of quantum non-locality many-body systems have are
essentially important issues. The study of quantum non-locality can complement
understanding the broad extensive study of quantum entanglement. Several
studies have been reported along this context, a quantum phase transition in
many-body system is characterized from the view of quantum non-locality in
some complex lattice models Campbell2010 ; Justino2012 ; Huang2013 ; Sun2014 ;
Lee2020 ; Wen2022 ; Liu2022  and disordered models Getelina2018 ; Liang2023 .
But there is no unified understanding. Also, there is an interesting
conjecture, two qubit quantum non-locality characterized by a violation of
CHSH inequality does not occur in (discrete) translational invariant systems
Oliveira2012 .

This work contributes to the understanding of the appearance of quantum non-
locality in many-body systems along the context discussed above. We study how
non-locality appears in the internal structure in a non-interacting spinless
fermion chain under quasi-periodic disorder, called extended Harper
modelHatsugai1990 ; Han1994 ; Liu2015 . For fermion many-particle states of
this model, two site and quantum non-locality is quantified by the CHSH
inequality and further multi-site quantum non-locality, which is three site
non-locality quantified by a Mermin–Klyshko-Svetlichny (MKS) polynomial
Svetlichny1987 ; Mermin1990 ; Belinskii ; Collins2002 ; Bancal2009 . We
clarify the internal structure of such quantum non-locality estimated by the
CHSH function and MKS polynomial, which has not been clarified so far.
Especially, due to $U(1)$ conservation (particle number conservation) the CHSH
and MKS polynomial can be efficiently calculated. We show how quantum non-
locality explicitly characterize non-interacting many-body fermion states and
its extended and critical phase transitions. In the globally averaged picture
of the maximum value of quantum non-locality is not broken, but within the
internal of the system the violation probabilities of the CHSH inequality and
MKS polynomial can become sufficiently finite.

The rest of this paper is organized as follows. In Sec. II, we introduce the
target fermion model and its basic properties. In Sec. III, the CHSH function
characterizing quantum non-locality between two sites and its practical
calculation method for this model are explained and numerical study is shown.
In Sec. IV, we further show the practical scheme of the calculation of MKS
polynomial applicable to this model and numerical study of adjacent three site
(qubit or spin) MKS polynomial. Section V is devoted for conclusion.

##  II Model

In this work, we consider the extended Harper model Hatsugai1990 ; Han1994 ,
the Hamiltonian of which is described by

$\displaystyle
H_{f}=\sum^{L-1}_{j=0}\biggr{[}J_{j}(f^{\dagger}_{j}f_{j+1}+\mbox{h.c.})+h_{j}n_{j}\biggl{]},$
(1)

where $f^{\dagger}_{j}$ and $f_{j}$ are fermion creation and annihilation
operators, the number operator is $n_{j}\equiv f^{\dagger}_{j}f_{j}$ and $L$
is the system size. Here, the hopping coupling and on-site potential are
quasi-periodic disorders,

$\displaystyle J_{j}$ $\displaystyle=$ $\displaystyle
1+\Delta_{1}\sin(2\pi\beta j+\phi_{0}),$ $\displaystyle h_{j}$
$\displaystyle=$ $\displaystyle\Delta_{2}\cos(2\pi\beta j+\phi_{0}),$

where $\beta=\frac{\sqrt{5}-1}{2}$ (golden ratio) and $\phi_{0}$ is a phase
shift. Throughout this work, we consider half filling (the total particle
number is $L/2$), open boundary condition and the couplings of the quasi-
periodic disorders are positive, $\Delta_{1(2)}\geq 0$.

So far, the phase diagram of the model of $H_{f}$ has been investigated in
Han1994 ; Liu2015 . The schematic phase diagram is shown in Fig. 1. There are
three phases: (i) Extended phase in the regime $\Delta_{1}\leq 1$ and
$\Delta_{2}\leq 2$, (ii) Critical phase in the regime $\Delta_{1}\geq 1$ and
$\Delta_{1}\geq\frac{1}{2}\Delta_{2}$, (iii) Localized phase in the regime
$\Delta_{2}\geq 2$ and $\Delta_{1}\leq\frac{1}{2}\Delta_{2}$. In the single
particle picture, these states are true for all single particle eigenstates.
There is no mobility edge in all phases Liu2015 . That is, if we consider
excited states of many-body case the many-body states are in the same phase
diagram. In particular, the critical phase includes broad correlation length
Han1994 . Its explicit characterization of this phase is a subtle but
important issue. In Ref.Liu2015 , a moderate value of participation ratio has
been employed, which qualitatively identify the phase. We shall show that
value of CHSH and MKS polynomials is a good indicator to identify the
transition point.

We study the quantum non-locality of the many-body fermion state with $L/2$
fermion particles conserved, the state of which is described by

$\displaystyle|\Psi\rangle=\prod^{L/2-1}_{\ell=0}\psi^{\dagger}_{i_{\ell}}|0\rangle,$
(2)

where $\psi^{\dagger}_{i_{\ell}}$ represents $i_{\ell}$-th single particle
eigenstate of $H_{f}$ with $i_{\ell}$-th lowest single particle energy,
described by
$\psi^{\dagger}_{i_{\ell}}=\sum^{L-1}_{j=0}c^{i_{\ell}}_{j}f^{\dagger}_{j}$
where $\sum^{L-1}_{j=0}|c^{i_{\ell}}_{j}|^{2}=1$ and $|0\rangle$ is a vacuum
state.

Figure 1: Schematic phase diagram of the extended Harper model. There are
three phases: extended, critical and localized phase. The phase boundary line
between critical and localized phases is $2\Delta_{1}=\Delta_{2}$.

##  III Investigation of quantum non-locality by CHSH function

Our interest is to investigate quantum non-local property in the non-
interecting many-body feremion state in the model of $H_{f}$. To this end, we
introduce the CHSH function, obtained from the correlations of two body spin
operators. In the non-interacting fermion, such spin operator can be
represented by Jordan-Wigner (JW) transformation LSM1961  shown in later.

###  III.1 CHSH function

The CHSH function is defined as

$\displaystyle
C(\vec{a}_{1},\vec{a}_{2},\vec{b}_{1},\vec{b}_{2})=\mathrm{tr}[\hat{\rho}_{\ell
m}\hat{M}_{2}].$ (3)

Here,

$\displaystyle\hat{M}_{2}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\biggr{[}\vec{a}_{1}\cdot\vec{\sigma}_{\ell}\otimes(\vec{b}_{1}+\vec{b}_{2})\cdot\vec{\sigma}_{m}$
(4)
$\displaystyle+\vec{a}_{2}\cdot\vec{\sigma}_{\ell}\otimes(\vec{b}_{1}-\vec{b}_{2})\cdot\vec{\sigma}_{m}\biggl{]},$

where $\ell,m$ represent a pair of spacial sites and $\vec{\sigma}_{\ell}$ is
a Pauli operator vector described by
$\vec{\sigma}_{\ell}=(\sigma^{x}_{\ell},\sigma^{y}_{\ell},\sigma^{z}_{\ell})$,
$\hat{\rho}_{\ell m}$ is a reduced density matrix on the two site pair
$(\ell,m)$ and the set of normalized three-dimensional vectors
$(\vec{a}_{1},\vec{a}_{2},\vec{b}_{1},\vec{b}_{2})$ acts as optimization
parameters. Then, $\hat{\rho}_{\ell m}$ is a reduced two site density matrix,
given by a spin expectation value and spin-spin correlations Horodecki1995 ,

$\displaystyle\hat{\rho}_{\ell m}$ $\displaystyle=$
$\displaystyle\frac{1}{4}\biggl{[}I\otimes
I+(\vec{r}\cdot\vec{\sigma}_{\ell})\otimes
I+I\otimes(\vec{s}\cdot\vec{\sigma}_{m})$ (5) $\displaystyle+$
$\displaystyle\sum_{\alpha,\beta=x,y,z}t_{\alpha\beta}\sigma^{\alpha}_{\ell}\otimes\sigma^{\beta}_{m}\biggr{]},$

where $I$ is a three dimensional identity matrix, $\vec{r}$, $\vec{s}$, and
$t_{\alpha\beta}$ can be obtained by their expectation values of the many-body
fermion state such as
$\vec{r}=(\langle\sigma^{x}_{\ell}\rangle,\langle\sigma^{y}_{\ell}\rangle,\langle\sigma^{z}_{\ell}\rangle)$,
$\vec{s}=(\langle\sigma^{x}_{m}\rangle,\langle\sigma^{y}_{m}\rangle,\langle\sigma^{z}_{m}\rangle)$
and $t_{\alpha\beta}=\langle\sigma^{\alpha}_{\ell}\sigma^{\beta}_{m}\rangle$
since $\mathrm{tr}_{(\ell,m)}[\hat{\rho}_{\ell m}\hat{X}_{\ell
m}]=\mathrm{tr}[\hat{\rho}_{0}\hat{X}_{\ell m}]$, in which $\hat{\rho}_{0}$ is
a density matrix constructed by many-body state Cheong2004  and $\hat{X}_{\ell
m}$ is an operator (observable) acting only on the sites, $\ell$ and $m$.

###  III.2 Non-local violation regime of CHSH function

As for a work Horodecki1995 , the maximum value of CHSH function given by Eq.
(3) denoted by $C_{max}$ can be obtained without optimization for the set of
the vectors $(\vec{a}_{1},\vec{a}_{2},\vec{b}_{1},\vec{b}_{2})$,

$\displaystyle C_{max}=2\sqrt{\lambda_{1}+\lambda_{2}},$ (6)

where $\lambda_{1}$ and $\lambda_{2}$ are the largest and second largest
eigenvalue of a $3\times 3$ matrix, $U=(T)^{t}T$ with
$(T)_{\alpha,\beta}=t_{\alpha\beta}$. From the spin picture of the model of
$H_{f}$, the spin-rotational symmetry in $x$-$y$ spin plane exists, then these
correlation functions $t_{\alpha\beta}$ satisfy
$\langle\sigma^{x}_{\ell}\sigma^{x}_{m}\rangle=\langle\sigma^{y}_{\ell}\sigma^{y}_{m}\rangle\neq\langle\sigma^{z}_{\ell}\sigma^{z}_{m}\rangle$
come1 , and also $t_{\alpha\beta}=0$ for $\alpha\neq\beta$, straightforwardly
obtained by JW fermion mapping and constraint of particle number conservation.
Then, the matrix $U$ is diagonal. Thus, the maximum value of the CHSH is given
by

$\displaystyle
C_{max}=\max\\{2\sqrt{2}|t_{xx}|,2\sqrt{|t_{zz}|^{2}+|t_{xx(yy)}|^{2}}\\}.$
(7)

Here, in quantum theory the maximal bound is known to be $C_{max}\leq
2\sqrt{2}$ Cirelson . This is consistent to that the maximal correlation value
$|t_{\alpha\beta}|\leq 1$. On the other hand, in the classical local theory (a
hidden variable classical theory) CHSH1969 , the maximum value of
$C(\vec{a}_{1},\vec{a}_{2},\vec{b}_{1},\vec{b}_{2})$ is bounded by
$\max|C(\vec{a}_{1},\vec{a}_{2},\vec{b}_{1},\vec{b}_{2})|\leq 2$. Thus, in
this work, we regard the regime $2<C_{max}\leq 2\sqrt{2}$ as a quantum non-
local violation regime. Our aim is to identify model-parameter regimes
appearing the non-local violation regime by the CHSH function in the model of
$H_{f}$.

###  III.3 Calculation method of correlation functions in non-interacting
many-body fermion states

The CHSH function can be obtained by spin-spin correlation functions.
According to the celebrated paper LSM1961 , the correlation functions can be
straightforwardly treated in fermion picture from a set of single particle
eigenstates, $\\{(\epsilon_{k},\vec{\psi}^{k})\\}$ with $k=0,1,\cdots,L-1$,
where $\vec{\psi}^{k}$ is $k$-th eigenstate of $H_{f}$, which is a real
$L$-component column vector and $\epsilon_{k}$ $k$-th eigenenergy. The spin-
spin correlation function can be obtained by employing the calculation method
proposed by LSM1961 ; Sachdev1999 ; Osborne2002 . Here we make use of the
particle conserving constraint.

At first, we consider the $z$-component spin-spin correlation for a pair of
site $(\ell,m)$. By using the JW transformation,

$\displaystyle\sigma^{x}_{j}$ $\displaystyle=$ $\displaystyle
f^{\dagger}_{j}e^{i\pi\sum^{j-1}_{k=0}n_{k}}+e^{-i\pi\sum^{j-1}_{k=0}n_{k}}f_{j},$
$\displaystyle\sigma^{y}_{j}$ $\displaystyle=$
$\displaystyle-i\biggr{[}f^{\dagger}_{j}e^{i\pi\sum^{j-1}_{k=0}n_{k}}+e^{-i\pi\sum^{j-1}_{k=0}n_{k}}f_{j}\biggl{]},$
$\displaystyle\sigma^{z}_{j}$ $\displaystyle=$ $\displaystyle 2n_{j}-1,$ (8)

the correlation is described by

$\displaystyle\langle\sigma^{z}_{\ell}\sigma^{z}_{m}\rangle$ $\displaystyle=$
$\displaystyle[G_{\ell,\ell}G_{m,m}-G_{m,\ell}G_{\ell,m}],$ (9)

where $G_{m,m}$ is the lattice Green function Cheong2004 ; Peschel2009 ,

$\displaystyle G_{i,j}\equiv\begin{cases}\langle
f^{\dagger}_{i}f_{j}\rangle+\langle f^{\dagger}_{j}f_{i}\rangle&\text{($i\neq
j$)}\\\ 2\langle f^{\dagger}_{i}f_{j}\rangle-1&\text{($i=j$)},\end{cases}$
(10)

and

$\displaystyle\langle f^{\dagger}_{i}f_{j}\rangle$ $\displaystyle=$
$\displaystyle\sum_{q\in\Lambda_{L/2}}\vec{\psi}^{q}_{i}(\vec{\psi}^{q}_{j})^{\dagger}.$
(11)

$\Lambda_{L/2}$ represent a set of $L/2$ single particle eigenstates. Note
that $G_{\ell,m}$ is symmetric, $G_{\ell,m}=G_{m,\ell}$ since the eigenvector
$\vec{\psi}^{k}$ is real for any $k$.

Other spin components of correlations are similarly given by

$\displaystyle\langle\sigma^{x}_{\ell}\sigma^{x}_{m}\rangle$ $\displaystyle=$
$\displaystyle\mathrm{det}\tilde{G}(\ell,m),$ (12)
$\displaystyle\tilde{G}(\ell,m)$ $\displaystyle\equiv$
$\displaystyle\begin{pmatrix}G_{\ell,\ell+1}&G_{\ell,\ell+2}&\dots&G_{\ell,m}\\\
G_{\ell+1,\ell+1}&G_{\ell+1,\ell+2}&\dots&G_{\ell+1,m}\\\
\vdots&\vdots&\ddots&\vdots\\\
G_{m-1,\ell+1}&G_{m-1,\ell+2}&\dots&G_{m-1,m}\end{pmatrix},$

and

$\displaystyle\langle\sigma^{y}_{\ell}\sigma^{y}_{m}\rangle$ $\displaystyle=$
$\displaystyle\mathrm{det}\bar{G}(\ell,m),$ (13)
$\displaystyle\bar{G}(\ell,m)$ $\displaystyle\equiv$
$\displaystyle\begin{pmatrix}G_{\ell+1,\ell}&G_{\ell+1,\ell+1}&\dots&G_{\ell,m-1}\\\
G_{\ell+2,\ell}&G_{\ell+2,\ell+1}&\dots&G_{\ell+1,m-1}\\\
\vdots&\vdots&\ddots&\vdots\\\
G_{m,\ell}&G_{m,\ell+1}&\dots&G_{m,m-1}\end{pmatrix}.$

Note that since the model of $H_{s}$ is not translational invariant, the above
spin-spin correlation functions are not described only by the distance of the
pair of sites $r=|\ell-m|$.

###  III.4 Numerical observation for CHSH function

We numerically observe the detailed behavior of the CHSH function, and treat
many different samples with different uniform random phase shift $\phi_{0}$.
Each sample is labeled by $s$ and the total number of samples is denoted by
$N_{s}$.

In this numerical study, the following quantities are considered. The first is
an averaged value of the maximum value of CHSH defined as

$\displaystyle\langle
C_{max}\rangle(r)\equiv\frac{1}{N_{s}}\sum_{s}\biggl{[}\frac{1}{L-r}\sum^{L-r-1}_{j=0}C^{s}_{max}(j,j+r)\biggr{]}.$
(14)

Here $C^{s}_{max}(j,j+r)$ is the value of CHSH function of Eq. (6) of a sample
$s$, where the pair of site is $j$ and $j+r$, $r$ is a distance between the
pair. From the above expression, the value of $\langle C_{max}\rangle(r)$
represents a global property of the degree of non-locality in the entire
system.

As the second quantity, we introduce the probability of the local CHSH
violation event denoted by $p_{v}$. This quantity is defined as

$\displaystyle p_{v}(r)$ $\displaystyle\equiv$
$\displaystyle\frac{1}{N_{s}}\sum_{s}v_{s}(r),$ (15) $\displaystyle v_{s}(r)$
$\displaystyle\equiv$
$\displaystyle\frac{1}{L-r}\sum^{L-r-1}_{j=0}H_{0}(C^{s}_{max}(j,j+r)-2).$
(16)

Here $H_{0}(x)$ is the Heaviside step function where $H_{0}(x)=1$ for $x>0$,
otherwise, $H_{0}(x)=0$. That is, if the CHSH function a pair of sites in the
system is in the non-local violation regime, the event is counted by
$H_{0}(x)$. The quantity $p_{v}(r)$ evaluates whether CHSH inequality breaking
is occurring somewhere in a two-site pair within the system. It gives
information on the CHSH inequality breaking from the local point view, giving
a different point of view from the globally average value of the CHSH function
of Eq. (14).

In the numerical calculation, two different many-body states are considered:
(i) $L/2$ particle ground state, corresponding to the set
$\\{i_{\ell}\\}=\\{i_{0},i_{1},\cdots,i_{L/2-1}\\}=\\{0,1,2,\cdots,L/2-1\\}$
in the state definition of Eq. (2), denoted as $|\psi_{g}\rangle$, (ii) an
typical excited state, set by
$\\{i_{\ell}\\}=\\{i_{0},i_{1},\cdots,i_{L/2-1}\\}=\\{L/4,L/4+1,L/4+2,\cdots,3L/4-1\\}$
in the state definition of Eq. (2), denoted as $|\psi_{e}\rangle$ and we focus
on two parameter sweeps as shown in Fig. 1: (A) varying $\Delta_{1}$ with
$\Delta_{2}=1$ corresponding to the parameter sweep between the extended and
critical phases, (B) varying $\Delta_{1}$ with the condition
$\Delta_{2}=2\Delta_{1}$, corresponding to the parameter sweep between the
extended and phase boundary line of the critical and localized phases.

Figure 2 shows the numerical results for the average value $\langle
C_{max}\rangle$ and the internal violation probability $p_{v}$ for the
groundstate $|\psi_{g}\rangle$ and the excited state $|\psi_{e}\rangle$, where
we set the parameter sweep (A) and observe the distance of two site pair up to
$r=6$.

Figure 2:  (a) Average value $\langle C_{max}\rangle$ for the groundstate
$|\psi_{g}\rangle$. (b) Internal violation probability $p_{v}$ for the
groundstate $|\psi_{g}\rangle$. (c) Average value $\langle C_{max}\rangle$ for
the excited state $|\psi_{e}\rangle$. (d) Internal violation probability
$p_{v}$ for the excited state $|\psi_{e}\rangle$. For all data, the parameter
sweep is of (A), where $\Delta_{2}=1$, and the distance of two site pair is
chosen by $r=1-6$. The system size is $L=64$ and the number of $\phi_{0}$
samples is $7\times 10^{3}$.

See Fig. 2(a), we observe that $\langle C_{max}\rangle$ exhibits a drop around
$\Delta_{1}=1$ for any $r$ for the state $|\psi_{g}\rangle$, indicating the
presence of the phase transition between the extended and critical phases,
consistent to the previous works Han1994 ; Liu2015 . That is, $\langle
C_{max}\rangle$ is a good indicator to identify a phase transition, and we
also find that the entire value of $\langle C_{max}\rangle$ on the parameter
sweep does not break the CHSH inequality, any non-local violation does not
appear globally in the system for any distance $r$, both phases and transition
points and $r=1$ value of $\langle C_{max}\rangle$ is larger than any other
$r$ case. We turn to observe $p_{v}$ as Fig. 2(b). Interestingly, a finite
value of $p_{v}$ ($\mathcal{O}(10^{-1})$) appears for $r=1$ case from the
regime near to critical phase but not for $r>1$. This implies that the CHSH
inequality breaks down with a finite probability in a portion of the internal
of the system. We also note that the violation probability $p_{v}$ in the
critical phase is larger than that of extended phase but $\langle
C_{max}\rangle$ in the critical phase is smaller than that of extended phase.
In addition, in Fig. 2(b), $p_{v}$ exhibits a much small finite value for
$r=3,5$ (generally, $r$ is in the odd case). The data implies the non-local
violation may occur as a rare event in the deep critical phase.

We next focus on the case of the excited state $|\psi_{e}\rangle$ for the
parameter sweep (A). In Fig. 2(c), $\langle C_{max}\rangle$ exhibits sudden
increase around $\Delta_{1}=1$ for any $r$, indicating the presence of the
phase transition of the excited state $|\psi_{e}\rangle$ between the extended
and critical phases. The entire value of $\langle C_{max}\rangle$ on the
parameter sweep indicates that the CHSH inequality does not break down
globally. However, we turn to observe $p_{v}$ as Fig. 2(d). Interestingly, we
find that a peak of a finite value of $p_{v}$ ($\mathcal{O}(10^{-1})$) appears
for $r=1$ and also for $r=2$ cases a finite value ($\mathcal{O}(10^{-2})$)
appears. This implies that the non-local quantum correlation for $r=1$ and $2$
occurs as a breakdown of the CHSH inequality with a finite probability in a
portion of the internal system around the phase transition. We also note that
the violation probability $p_{v}$ for $r=2$ once exhibit a sudden increase
around $\Delta_{1}=1$ and then decreases suddenly after passing $\Delta_{1}=1$
and increase again with increase $\Delta_{1}$. This peculiar behavior suggests
that the phase change between the extended and critical phases may not be a
single phase transition. A similar behavior appears in the later calculation
of the MKS polynomial.

Figure 3: (a) Average value $\langle C_{max}\rangle$ for the groundstate
$|\psi_{g}\rangle$. (b) Internal violation probability $p_{v}$ for the
groundstate $|\psi_{g}\rangle$. (c) Average value $\langle C_{max}\rangle$ for
the excited state $|\psi_{e}\rangle$. (d) Internal violation probability
$p_{v}$ for the excited state $|\psi_{e}\rangle$. For all data, the parameter
sweep is of (B) and the distance of two site pair is chosen by $r=1-6$. The
system size is $L=64$ and the number of $\phi_{0}$ samples is $7\times
10^{3}$.

We move to the case of the parameter sweep (B). The behavior of $\langle
C_{max}\rangle$ for the groundstate $|\psi_{g}\rangle$ is displayed in Fig.
3(a). All data exhibits a drop around $\Delta_{1}=1$ for any $r$ for the state
$|\psi_{g}\rangle$, indicating the presence of the phase change of the ground
state between the extended phase and the regime on the phase boundary between
the critical and localized phases, again consistent to the previous works
Han1994 ; Liu2015 . $\langle C_{max}\rangle$ well characterize the vanishing
extended phase. The entire value of $\langle C_{max}\rangle$ on the parameter
sweep indicates that the CHSH inequality does not break down globally in the
system for any distance $r$. Also, the value of $\langle C_{max}\rangle$ for
$r\geq 2$ is same on the boundary, implying a strong criticality, that is, the
degree of correlation exhibits scale free tendency. We further observe $p_{v}$
for the groundstate $|\psi_{g}\rangle$ as Fig. 3(b). Interestingly, a large
finite value of $p_{v}$ ($\mathcal{O}(10^{-1})$) appears only for $r=1$ case
but not for $r>1$ and the regime starts around a regime vanishing the extended
phase $\Delta_{1}\sim 1$. Only the nearest neighbor correlation exhibits
quantum non-locality as the break down of the CHSH inequality with a finite
probability in a portion of the internal of the system on the phase boundary
between the critical and localized phases.

We further observe the case of the excited state $|\psi_{e}\rangle$. In Fig.
3(c), $\langle C_{max}\rangle$ exhibits a sudden increase around
$\Delta_{1}=1$ for any $r$, indicating the presence of vanishing the extended
phase in the excited state $|\psi_{e}\rangle$. The entire value of $\langle
C_{max}\rangle$ on the parameter sweep indicates that the CHSH inequality does
not break down globally in the system for any distance $r$. However, we find
that as shown in Fig. 3(d), a finite value of $p_{v}$ ($\mathcal{O}(10^{-1})$)
appears for $r=1$ and $r=2$ cases after stating to vanish the extended phase.
This implies that the correlation for $r=1$ and $2$ exhibits quantum
nonlocality as the break down of the CHSH inequality with a finite probability
in a portion of the internal system on the phase boundary between the critical
and localized phases.

We also care the system size dependence of $\langle C_{max}\rangle$ and
$p_{v}$. The properties of the transition and its entire behavior are almost
same. The results are shown in Appendix C.

##  IV Investigation of three-qubit MKS polynomial

We next observe an extension form of the CHSH function. The general form is
called Mermin–Klyshko-Svetlichny (MKS) polynomial Svetlichny1987 ; Mermin1990
; Belinskii ; Belinskii ; Collins2002 , which can characterize $n$-qubit
multipartite quantum non-locality, also can capture its hierarchical structure
of non-locality Wen2022 . The brief explanation of general form of the MKS
polynomial is in Appendix A.

In this numerical study, we limit adjacent three-qubit MKS (3MKS) polynomial.
The concrete form for a three-qubit subsystem (three spatial sites) we
consider here is given as the expectation value of 3MKS polynomial
Svetlichny1987 ; Mermin1990 ; Belinskii ; Collins2002

$\displaystyle C^{3MKS}({\bf
a})=\frac{1}{2}\biggr{[}\mathrm{tr}_{3}(\hat{\rho}_{3}\hat{M}_{3})+\mathrm{tr}_{3}(\hat{\rho}_{3}\hat{M}^{\prime}_{3})\biggl{]},$
(17)

where ${\bf
a}=(\vec{a}_{j_{1}},\vec{a}^{\prime}_{j_{1}},\vec{a}_{j_{2}},\vec{a}^{\prime}_{j_{2}},\vec{a}_{j_{3}},\vec{a}^{\prime}_{j_{3}})$,
corresponds to the set of optimization vector parameters, $\hat{\rho}_{3}$ is
an adjacent three-site reduced density matrix and $\hat{M}^{(^{\prime})}_{3}$
is 3MKS operator, constituted by some lower-lank operators as

$\displaystyle\hat{M}_{3}=\frac{1}{2}\hat{M}_{2}\otimes(\vec{a}_{j_{3}}+\vec{a}^{\prime}_{j_{3}})\cdot\vec{\sigma}_{j_{3}}+\frac{1}{2}\hat{M}^{\prime}_{2}\otimes(\vec{a}_{j_{3}}-\vec{a}^{\prime}_{j_{3}})\cdot\vec{\sigma}_{j_{3}},$
$\displaystyle\hat{M}^{\prime}_{3}=\frac{1}{2}\hat{M}^{\prime}_{2}\otimes(\vec{a}_{j_{3}}+\vec{a}^{\prime}_{j_{3}})\cdot\vec{\sigma}_{j_{3}}+\frac{1}{2}\hat{M}_{2}\otimes(\vec{a}^{\prime}_{j_{3}}-\vec{a}_{j_{3}})\cdot\vec{\sigma}_{j_{3}},$
$\displaystyle\hat{M}_{2}=\frac{1}{2}\hat{M}_{1}\otimes(\vec{a}_{j_{2}}+\vec{a}^{\prime}_{j_{2}})\cdot\vec{\sigma}_{j_{2}}+\frac{1}{2}\hat{M}^{\prime}_{1}\otimes(\vec{a}_{j_{2}}-\vec{a}^{\prime}_{j_{2}})\cdot\vec{\sigma}_{j_{2}},$
$\displaystyle\hat{M}^{\prime}_{2}=\frac{1}{2}\hat{M}^{\prime}_{1}\otimes(\vec{a}_{j_{2}}+\vec{a}^{\prime}_{j_{2}})\cdot\vec{\sigma}_{j_{2}}+\frac{1}{2}\hat{M}_{1}\otimes(\vec{a}^{\prime}_{j_{2}}-\vec{a}_{j_{2}})\cdot\vec{\sigma}_{j_{2}},$
$\displaystyle\hat{M}_{1}=\vec{a}_{j_{1}}\cdot\vec{\sigma}_{j_{1}},\>\>\hat{M}^{\prime}_{1}=\vec{a}^{\prime}_{j_{1}}\cdot\vec{\sigma}_{j_{1}},$

where each optimization vector $\vec{a}^{(^{\prime})}_{j}$ is a spherical unit
vector parameterized by two angle parameters as
$\vec{a}^{(^{\prime})}_{j}=(\cos\theta^{(^{\prime})}_{j}\sin\phi^{(^{\prime})}_{j},\sin\theta^{(^{\prime})}_{j}\sin\phi^{(^{\prime})}_{j},\cos\phi^{(^{\prime})}_{j})$
and the label $(j_{1},j_{2},j_{3})$ represents a three adjacent lattice sites
in the system. Our interest is to estimate a maximum value of $C^{3MKS}$
denoted by $C^{3MKS}_{max}$. To this end, we solve the optimization problem
for 12 rotational angle parameters by employing a numerical Scipy package
Scipy_opt .

In the local classical theory (a local hidden variable theory), the maximum
value of $C^{3MKS}_{max}$ is bounded by $C^{3MKS}_{max}\leq 1$ Gisin1998 ;
Collins2002 . On the other side, in quantum theory, the maximum value of
$C^{3MKS}_{max}$ is bounded by $C^{3MKS}_{max}\leq\sqrt{2}$, where
$C^{3MKS}_{max}=\sqrt{2}$ indicates full three-body non-locality among the
three site. Thus, the quantum non-locality for the three-sites is
characterized by $1<C^{3MKS}_{max}\leq\sqrt{2}$, called 3MKS non-local
violation regime.

Figure 4: (a) Average value $\langle C^{3MKS}_{max}\rangle$ for the parameter
sweep (A). (b) Internal violation probability $p^{MKS}_{v}$ for the parameter
sweep (A). (c) Average value $\langle C^{3MKS}_{max}\rangle$ for the parameter
sweep (B). (d) Internal violation probability $p^{MKS}_{v}$ for the parameter
sweep (B). The system size is $L=64$ and the number of $\phi_{0}$ samples is
$8\times 10^{3}$. The labels “GS” and “EX” represent the states
$|\psi_{g}\rangle$ and $|\psi_{e}\rangle$.

As an example, in only three qubit system, the GHZ state $\displaystyle{|{\rm
GHZ}\rangle=\frac{1}{\sqrt{2}}[|111\rangle+|000\rangle]}$ is known to give a
maximum quantum non-locality, $C^{3MKS}_{max}=\sqrt{2}$ Collins2002 .

Let us move to numerical observations about the detailed behavior of the 3MKS,
as same to the case of the CHSH function. We set the three adjacent sites
around the center of the system, $j_{1}=L/2-1$, $j_{2}=L/2$, and
$j_{3}=L/2+1$, and treat many different samples with different $\phi_{0}$. We
explain the detail of the calculation of 3MKS for the fermion many-body system
as shown in Appendix B. We observe the averaged value of the maximum value of
3MKS,

$\displaystyle\langle
C^{3MKS}_{max}\rangle\equiv\frac{1}{N_{s}}\sum_{s}C^{3MKS,s}_{max},$ (18)

where we denote 3MKS of a sample $s$ by $C^{3MKS,s}_{max}$, and also observe
the violation event probability

$\displaystyle p^{MKS}_{v}$ $\displaystyle\equiv$
$\displaystyle\frac{1}{N_{s}}\sum_{s}v^{MKS}_{s},$ (19)

with $v^{MKS}_{s}\equiv H_{0}(C^{3MKS,s}_{max}-1)$. As same to the case of the
calculation of the CHSH function, we pick up the two parameter sweeps (A) and
(B), and focus on the two many-body states, $|\psi_{g}\rangle$ and
$|\psi_{e}\rangle$.

The behavior of $\langle C^{3MKS}_{max}\rangle$ for the parameter sweep (A) is
displayed in Fig. 4(a). The data of $|\psi_{g}\rangle$ exhibits a peak at
$\Delta_{1}=1$ and the data of $|\psi_{e}\rangle$ exhibits a drop at
$\Delta_{1}=1$, indicating the presence of the phase transition in the
different behavior. $\langle C^{3MKS}_{max}\rangle$ well characterize the
phase transition. The entire value of $\langle C^{3MKS}_{max}\rangle$ on the
parameter sweep does not enter into the 3MKS non-local violation regime. But,
we further observe $p^{3MKS}_{v}$ for the parameter sweep (A) as shown in Fig.
4(b). We find that for the state $|\psi_{g}\rangle$, a sudden increase of
$p^{3MKS}_{v}$ up to $\mathcal{O}(10^{-1})$ appears around $\Delta_{1}=1$ and
the 3MKS exhibits the non-local violation regime with a finite probability in
the internal of the system. We move to observe the data for the state
$|\psi_{e}\rangle$. We find that for the state $|\psi_{e}\rangle$,
$p^{3MKS}_{v}$ starts to have a finite value up to $\mathcal{O}(10^{-1})$ from
entering the critical phase and the 3MKS exhibits the non-local violation
regime with a finite probability in the internal of the system. Further, note
that a large fluctuation behavior around $\Delta_{1}=1$ is observed. This
behavior is similar to the behavior of the CHSH function shown in $r=2$ data
of Fig. 2(d), indicating the phase change between the extended and critical
phases may not be a single phase transition.

Finally, we turn to the results for the parameter sweep (B). In Fig. 4(c),
$\langle C_{max}\rangle$ exhibits similar behavior to the case of the
parameter sweep (A) and characterize vanishing the extended phase for both
states $|\psi_{g}\rangle$ and $|\psi_{e}\rangle$. The entire value of $\langle
C_{max}\rangle$ on the parameter sweep does not exhibit the MKS non-local
violation regime. As shown in Fig. 4(d), the event of the MKS non-local
violation occurs only around the phase change as observed that only for the
state $|\psi_{g}\rangle$, a finite value of $p_{v}$ ($\mathcal{O}(10^{-2})$)
appears at $\Delta_{1}=1$. This indicate that the event of the MKS non-local
violation occurs around the tricritical point of the model in the ground state
$|\psi_{g}\rangle$ with a finite probability.

##  V Conclusion

We showed that the CHSH function and MKS polynomial are good indicators to
characterize the phase transitions in disordered fermion systems. As an
example, we studied non-local properties in the extended Harper model. We
found that the non-local violation regime of both the CHSH function and 3MKS
locally appears in a critical phase and around its transition points with a
finite probability. Such a violation does not appear as an averaged value.
Further, the CHSH function and 3MKS exhibits similar behavior and they are
also useful for detecting a phase transition for excited states. These
quantities can be useful for understanding the detailed structure of some
complex phase transitions in disordered systems.

We expect that this characterization approach suggested in this work has broad
applications for various non-interacting systems, such as long-range hopping
model, e.g. how the long-range nature changes properties of quantum non-
locality in many-body systems.

## ACKNOWLEDGMENTS

The work is supported by JSPS KAKEN-HI Grant Number 23K13026 (Y.K.), Japan.

## Appendix

## A. General MKS polynomial

General MKS polynomials for $n$-qubits can be recursively defined depending on
even or odd $n$ Collins2002 . The two types of MKS polynomias are introduced
as the following forms

$\displaystyle S^{MKS,1}_{n}({\bf a})$ $\displaystyle=$
$\displaystyle\mathrm{tr}(\hat{\rho}_{n}\hat{M}_{n}),$ (A.1 ) $\displaystyle
S^{MKS,2}_{n}({\bf a})$ $\displaystyle=$
$\displaystyle\mathrm{tr}(\hat{\rho}_{n}\hat{M}^{{}^{\prime}}_{n}).$ (A.2 )

Here $\hat{\rho}_{n}$ is a reduced density matrix for $n$-qubit subsystem in a
system, ${\bf
a}=(\vec{a}_{j_{1}},\vec{a}^{\prime}_{j_{1}},\vec{a}_{j_{2}},\vec{a}^{\prime}_{j_{2}},\cdots,\vec{a}_{j_{n}},\vec{a}^{\prime}_{j_{n}})$
including total $4n$ optimization parameters, and the above $n$-qubit MKS
operator $\hat{M}_{n}$ is recursively defined by Gisin1998

$\displaystyle\hat{M}_{n}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\hat{M}_{n-1}\otimes(\vec{a}_{j_{n}}+\vec{a}^{{}^{\prime}}_{j_{n}})\cdot{\vec{\sigma}}_{j_{n}}$
(A.3 )
$\displaystyle+\frac{1}{2}\hat{M}^{{}^{\prime}}_{n-1}\otimes(\vec{a}_{j_{n}}-\vec{a}^{{}^{\prime}}_{j_{n}})\cdot{\vec{\sigma}}_{j_{n}},$

where the operator $\hat{M}^{{}^{\prime}}_{n}$ is obtained from $\hat{M}_{n}$
by exchanging all the primed and non-primed $\vec{a}$’s and $n=1$ case are
given by $\hat{M}_{1}=\vec{a}_{j_{1}}\cdot\vec{\sigma}_{j_{1}}$ and
$\hat{M}^{\prime}_{1}=\vec{a}^{\prime}_{j_{1}}\cdot\vec{\sigma}_{j_{1}}$.

As for Ref.Collins2002 , as a measure of quantum non-locality, the form of the
generalized MKS polynomial depends on whether $n$ is even or odd as

$\displaystyle C^{nMKS}\equiv\begin{cases}S^{MKS,1}_{n}({\bf
a})&\text{($n=$even)}\\\ \displaystyle{\frac{1}{2}\biggr{[}S^{MKS,1}_{n}({\bf
a})+S^{MKS,2}_{n}({\bf a})\biggl{]}}&\text{($n=$odd)}\end{cases}.$ (A.4 )

The maximum optimized value denoted by $C^{nMKS}_{max}$ ($\geq 0$) of the
above can evaluates $n$-qubit (or $n$-site) non-locality. If the theory is
local and classical (a local hidden variable model), $C^{nMKS}_{max}\leq 1$.
On the other hand, if the theory exhibits a quantum non-locality, the above
quantity $C^{nMKS}_{max}$ is in regimes; $1<C^{nMKS}_{max}\leq\sqrt{2^{n-1}}$
for even $n$ and $1<C^{nMKS}_{max}\leq\sqrt{2^{n-2}}$ for odd $n$. That is,
the maximum bound values for the quantum non-locality are $\sqrt{2^{n-1}}$ for
even $n$ and $\sqrt{2^{n-2}}$, respectively.

In this work, for $n=2$ case, we observe $C_{max}\equiv 2C^{nMKS}_{max}$ by
multiplying a factor $2$ respecting the form of CHSH inequality CHSH1969 ,
while we use the form of $C^{nMKS}_{max}$ of the MKS polynomial for $n=3$
case.

## B. Detail of calculation of $n=3$ MKS polynomial from fermion picture

The 3MKS of Eq. (17) is represented by a sum of three-spin correlation
functions,

$\displaystyle C^{3MKS}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\biggl{[}\mathrm{tr}(\hat{\rho}_{3}\hat{M}_{2}\otimes(\vec{a}^{{}^{\prime}}_{j_{3}}\cdot{\vec{\sigma}}_{j_{3}}))$
$\displaystyle+\mathrm{tr}(\hat{\rho}_{3}\hat{M}^{\prime}_{2}\otimes(\vec{a}_{j_{3}}\cdot\vec{\sigma}_{j_{3}}))\biggr{]}$
$\displaystyle=$
$\displaystyle\frac{1}{2}\biggl{[}\mathrm{tr}(\hat{\rho}_{3}\hat{M}_{2}\otimes(\vec{a}^{{}^{\prime}}_{j_{3}}\cdot{\vec{\sigma}}_{j_{3}}))$
$\displaystyle+\mathrm{tr}(\hat{\rho}_{3}\hat{M}^{\prime}_{2}\otimes(\vec{a}_{j_{3}}\cdot\vec{\sigma}_{j_{3}}))\biggr{]}$
$\displaystyle=$
$\displaystyle\sum_{(\alpha,\beta,\gamma)\in\mathcal{S}}\frac{1}{2}\biggr{[}a^{\alpha}_{j_{1}}a^{\beta}_{j_{2}}a^{\prime\gamma}_{j_{3}}+a^{\prime\alpha}_{j_{1}}a^{\beta}_{j_{2}}a^{\prime\gamma}_{j_{3}}$
$\displaystyle+a^{\alpha}_{j_{1}}a^{\prime\beta}_{j_{2}}a^{\prime\gamma}_{j_{3}}-a^{\prime\alpha}_{j_{1}}a^{\prime\beta}_{j_{2}}a^{\prime\gamma}_{j_{3}}\biggl{]}\langle\sigma^{\alpha}_{j_{1}}\sigma^{\beta}_{j_{2}}\sigma^{\gamma}_{j_{3}}\rangle$
$\displaystyle+\sum_{(\alpha,\beta,\gamma)\in\mathcal{S}}\frac{1}{2}\biggr{[}a^{\prime\alpha}_{j_{1}}a^{\prime\beta}_{j_{2}}a^{\gamma}_{j_{3}}+a^{\alpha}_{j_{1}}a^{\prime\beta}_{j_{2}}a^{\gamma}_{j_{3}}$
$\displaystyle+a^{\prime\alpha}_{j_{1}}a^{\beta}_{j_{2}}a^{\gamma}_{j_{3}}-a^{\alpha}_{j_{1}}a^{\beta}_{j_{2}}a^{\gamma}_{j_{3}}\biggl{]}\langle\sigma^{\alpha}_{j_{1}}\sigma^{\beta}_{j_{2}}\sigma^{\gamma}_{j_{3}}\rangle,$

where $\langle\cdot\rangle$ means expectation value for the many-body state
$|\Psi\rangle$ of Eq. (2), the optimization parameters are

$\displaystyle\vec{a}_{j}(\theta_{j},\phi_{j})$ $\displaystyle\equiv$
$\displaystyle(a^{x}_{j},a^{y}_{j},a^{z}_{j})$ $\displaystyle=$
$\displaystyle(\cos\theta_{j}\sin\phi_{j},\sin\theta_{j}\sin\phi_{j},\cos\phi_{j}),$
$\displaystyle\vec{a}^{\prime}_{j}(\theta^{\prime}_{j},\phi^{\prime}_{j})$
$\displaystyle\equiv$ $\displaystyle(a^{\prime x}_{j},a^{\prime
y}_{j},a^{\prime z}_{j})$ $\displaystyle=$
$\displaystyle(\cos\theta^{\prime}_{j}\sin\phi^{\prime}_{j},\sin\theta^{\prime}_{j}\sin\phi^{\prime}_{j},\cos\phi^{\prime}_{j}),$

and

$\displaystyle\mathcal{S}=$ $\displaystyle\\{(z,z,z),(x,z,x),(z,x,x),$
$\displaystyle(x,x,z),(y,z,y),(y,y,z),(z,y,y)\\}.$

The above sums over the labels of spin component, $(\alpha,\beta,\gamma)$. The
restricted choice of $\mathcal{S}$ originates from the constraint of particle
number conservation in the fermion model.

Practically, each expectation values of
$\langle\sigma^{\alpha}_{j_{1}}\sigma^{\beta}_{j_{2}}\sigma^{\gamma}_{j_{3}}\rangle$
with $(\alpha,\beta,\gamma)\in\mathcal{S}$ can be calculated by transforming
them into the fermion picture though the Jordan-Winger transformation. The
concrete forms for correlations given by the element of $\mathcal{S}$ are
given as follows:

$\displaystyle\langle\sigma^{x}_{j_{1}}\sigma^{z}_{j_{2}}\sigma^{x}_{j_{3}}\rangle$
$\displaystyle=$ $\displaystyle\langle
f^{\dagger}_{j_{1}}e^{i\sum^{j_{1}-1}_{m=0}n_{m}}(2n_{j_{2}}-1)e^{-i\sum^{j_{3}-1}_{m=0}n_{m}}f_{j_{3}}\rangle$
$\displaystyle+$
$\displaystyle\langle(f_{j_{1}}e^{-i\sum^{j_{1}-1}_{m=0}n_{m}}(2n_{j_{2}}-1)e^{i\sum^{j_{3}-1}_{m=0}n_{m}}f^{\dagger}_{j_{3}})\rangle,$
$\displaystyle\langle\sigma^{x}_{j_{1}}\sigma^{x}_{j_{2}}\sigma^{z}_{j_{3}}\rangle$
$\displaystyle=$ $\displaystyle\langle
f^{\dagger}_{j_{1}}f_{j_{2}}(2n_{j_{3}}-1)\rangle-\langle
f_{j_{1}}f^{\dagger}_{j_{2}}(2n_{j_{3}}-1)\rangle,$
$\displaystyle\langle\sigma^{z}_{j_{1}}\sigma^{x}_{j_{2}}\sigma^{x}_{j_{3}}\rangle$
$\displaystyle=$
$\displaystyle\langle(2n_{j_{1}}-1)f^{\dagger}_{j_{2}}f_{j_{3}}\rangle-\langle(2n_{j_{1}}-1)f_{j_{2}}f^{\dagger}_{j_{3}}\rangle,$
$\displaystyle\langle\sigma^{z}_{j_{1}}\sigma^{z}_{j_{2}}\sigma^{z}_{j_{3}}\rangle$
$\displaystyle=$ $\displaystyle 8\langle n_{j_{1}}n_{j_{2}}n_{j_{3}}\rangle$
$\displaystyle-4\langle n_{j_{1}}n_{j_{2}}\rangle-4\langle
n_{j_{2}}n_{j_{3}}\rangle-4\langle n_{j_{3}}n_{j_{1}}\rangle$
$\displaystyle+2\langle n_{j_{1}}\rangle+2\langle n_{j_{2}}\rangle+2\langle
n_{j_{3}}\rangle-1.$

By using Wick’s theorem, each correlations can be further written as a
multiple form of the fermion correlation functions of a non-interecting many-
body state, $\langle f^{\dagger}_{i}f_{j}\rangle$ obtained by Eq. (11) in the
main text. Also, since the Hamiltonian of Eq. (1) has a symmetry of
replacement $x\longleftrightarrow y$ for the spin components, some
relationship like
$\langle\sigma^{x}_{j_{1}}\sigma^{z}_{j_{2}}\sigma^{x}_{j_{3}}\rangle=\langle\sigma^{y}_{j_{1}}\sigma^{z}_{j_{2}}\sigma^{y}_{j_{3}}\rangle$
exists. Also, this fact is straightforwardly checked only by representing them
in fermion picture under particle conservation.

Figure 5: System size dependence for average value $\langle C_{max}\rangle$
for $r=1$ [(a)] and internal violation probability $p_{v}$ for $r=1$ [(b)] in
the parameter sweep (A). System size dependence for average value $\langle
C_{max}\rangle$ for $r=1$ [(c)] and internal violation probability $p_{v}$ for
$r=1$ in the parameter sweep (B). For all data, the number of $\phi_{0}$
samples is $3\times 10^{3}$-$1\times 10^{4}$. Figure 6: System size dependence
for average value $\langle C^{3MKS}_{max}\rangle$ [(a)] and internal violation
probability $p^{MKS}_{v}$[(b)] in the parameter sweep (A). For all data, the
number of $\phi_{0}$ samples is $7\times 10^{3}$.

As the result, if we consider the adjacent three site case
$j_{1}=j_{2}-1=j_{3}-2$, then in this system we only need to calculate total
seven patterns of the three body spin correlations coming from the set of
$\mathcal{S}$ to calculate the 3MKS.

## C. System size dependence

In this appendix, we show the system size dependence of our numerics for some
typical cases.

Figures 5 (a) and (b) are the system size dependence of $\langle
C_{max}\rangle$ for $r=1$ and $p_{v}$ for $r=1$, where we observe the
groundstate $|\psi_{g}\rangle$ with the parameter sweep (A). We verified that
the dependence is entirely much small for $\langle C_{max}\rangle$ and $p_{v}$
exhibits tiny system size dependence in the extended phase, but no system size
dependence in the critical phase.

Figures 5 (c) and (d) are the system size dependence of $\langle
C_{max}\rangle$ for $r=1$ and $p_{v}$ for $r=1$, where we observe the
groundstate $|\psi_{g}\rangle$ with the parameter sweep (B). We verified that
the dependence is entirely much small for $\langle C_{max}\rangle$ and $p_{v}$
also exhibits tiny system size dependence in the extended phase, but no system
size dependence on the phase boundary.

Furthermore, figures 6 (a) and (d) are the system size dependence of $\langle
C^{3MKS}_{max}\rangle$, where we observe the groundstate $|\psi_{g}\rangle$
with the parameter sweep (A). We verified that the dependence of $\langle
C^{3MKS}_{max}\rangle$ is much small around the phase transition point and in
the critical phase, while in the extended phase the small system size
dependence appears. $p^{MKS}_{v}$ also exhibits a small system size dependence
in the critical phase, but no system size dependence around phase transition
point.

## References

  * (1) R. Horodecki, P. Horodecki, M. Horodecki and K. Horodecki, Rev. Mod. Phys. 81, 865 (2009).
  * (2) J. Eisert, M. Cramer, and M. B. Plenio, Rev. Mod. Phys. 82, 277 (2010).
  * (3) B. Zeng, X. Chen, D.-L. Zhou, and X.-G. Wen, Quantum Information Meets Quantum Matter: From Quantum Entanglement to Topological Phases of Many-Body Systems (Springer, 2019).
  * (4) J.S. Bell, Physics 1, 195 (1964).
  * (5) J. F. Clauser, M. A. Horne, A. Shimony, R. A. Holt, Phys. Rev. Lett. 23, 880 (1969).
  * (6) R. F. Werner, Phys. Rev. A 40, 4277 (1989).
  * (7) R. Augusiak, M. Demianowicz, J. Tura, A. Acín, Phys. Rev. Lett. 115, 030404 (2015).
  * (8) N. Gisin, Phys. Lett. A 154, 201 (1991)
  * (9) S. Campbell and M. Paternostro, Phys. Rev. A 82, 042324 (2010).
  * (10) L. Justino and T. R. De Oliveira, Phys. Rev. A 85, 052128 (2012).
  * (11) H. L. Huang, Z. Y. Sun, and B. Wang, Eur. Phys. J. B 86, 279 (2013).
  * (12) Z. Y. Sun, S. Liu, H. L. Huang, D. Zhang, Y. Y. Wu, J. Xu, B. F. Zhan, H. G. Cheng, C. B. Duan, and B. Wang, Phys. Rev. A 90, 062129 (2014).
  * (13) H. X. Wen, Z. Y. Sun, H. G. Cheng, D. Zhang, and Y. Y. Wu, Eur. Phys. J. B 95, 95:148 (2022).
  * (14) Y. Liu, M. Li, J. Bao, B. Guo, and Z. Sun, Phys. Lett. A 450, 128396 (2022).
  * (15) D. Lee and W. Son, Entropy 22, 1282 (2020).
  * (16) J. C. Getelina, T. R. de Oliveira, and J. A. Hoyos, Phys. Lett. A 382, 2799 (2018).
  * (17) Z. Liang, J. Bao, L. Shen, B. Guo, and Z. Sun, Phys. Lett. Sect. A 472, 128810 (2023).
  * (18) T. R. de Oliveira, A. Saguia, and M. S. Sarandy, EPL, 100, 60004 (2012).
  * (19) Y. Hatsugai and M. Kohmoto, Phys. Rev. B 42, 8282 (1990).
  * (20) J. H. Han, D. J. Thouless, H. Hiramoto, and M. Kohmoto, Phys. Rev. B 50, 11365 (1994).
  * (21) F. Liu, S. Ghosh, and Y. D. Chong, Phys. Rev. B 91, 014108 (2015).
  * (22) G. Svetlichny, Phys. Rev. D 35 3066 (1987).
  * (23) N. D. Mermin, Phys. Rev. Lett. 65, 1838 (1990)
  * (24) A. V. Belinskii, D. N. Klyshko, Physics-Uspekhi 36, 653 (1993).
  * (25) D. Collins, N. Gisin, S. Popescu, D. Roberts, and V. Scarani, Phys. Rev. Lett. 88, 1704051 (2002).
  * (26) J. D. Bancal, C. Branciard, N. Gisin, and S. Pironio, Phys. Rev. Lett. 103, 090503 (2009).
  * (27) E. Lieb, T. Schultz and D. Mattis, Ann. Phys. (New York) 16, 407 (1961).
  * (28) R. Horodecki, P. Horodecki, and M. Horodecki, Phys. Lett. A 200, 340 (1995).
  * (29) S. A. Cheong and C. L. Henley, Phys. Rev. B 69, 075111 (2004).
  * (30) A similar consideration has been used in Refs. Justino2012 ; Getelina2018 . In addition, we also verified the condition by using the numerical exact diagonalization for the model of $H_{f}$.
  * (31) B. S. Cirel’son, Lett. Math. Phys. 4, 93 (1980).
  * (32) S. Sachdev, Quantum Phase Transitions (Cambridge University Press, Cambridge, UK, 1999).
  * (33) T. J. Osborne and M. A. Nielsen, Phys. Rev. A 66, 032110 (2002).
  * (34) I. Peschel and V. Eisler, J. Phys. A: Math. Theor. 42, 504003 (2009).
  * (35) https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.minimize.html
  * (36) N. Gisin and H. Bechmann-Pasquinucci, Phys. Lett. A 246, 1 (1998).


# Do not Mask Randomly:
Effective Domain-adaptive Pre-training by Masking In-domain Keywords

Shahriar Golchin†, Mihai Surdeanu†, Nazgol Tavabi⋄, Ata Kiapour⋄
†University of Arizona, Tucson, AZ, USA
⋄Harvard Medical School, Boston, MA, USA
<EMAIL_ADDRESS>

###### Abstract

We propose a novel task-agnostic in-domain pre-training method that sits
between generic pre-training and fine-tuning. Our approach selectively masks
in-domain keywords, i.e., words that provide a compact representation of the
target domain. We identify such keywords using KeyBERT (Grootendorst, 2020).
We evaluate our approach using six different settings: three datasets combined
with two distinct pre-trained language models (PLMs). Our results reveal that
the fine-tuned PLMs adapted using our in-domain pre-training strategy
outperform PLMs that used in-domain pre-training with random masking as well
as those that followed the common pre-train-then-fine-tune paradigm. Further,
the overhead of identifying in-domain keywords is reasonable, e.g., 7–15% of
the pre-training time (for two epochs) for BERT Large (Devlin et al.,
2019).111The code for all of our experiments is available at
https://github.com/shahriargolchin/do-not-mask-randomly.

##  1 Introduction

Employing large pre-trained language models (PLMs) is currently a common
practice for most natural language processing (NLP) tasks (Tunstall et al.,
2022). A two-stage pre-train-then-fine-tune framework is usually used to
adapt/fine-tune PLMs to downstream tasks (Devlin et al., 2019). However,
motivated by ULMFiT (Howard and Ruder, 2018) and ELMo (Peters et al., 2018),
Gururangan et al. (2020) showed that incorporating in-domain pre-training
(also known as domain-adaptive pre-training) between generic pre-training and
fine-tuning stages can lead to further performance improvements in downstream
tasks because it “pulls” the PLM towards the target domain. At this
intermediate stage, the domain adaptation for PLMs is typically handled by
continuing pre-training in the same way, i.e., using randomly-masked tokens on
unstructured in-domain data Devlin et al. (2019). Here, we argue that this
intermediate pre-training should be performed differently, i.e., masking
should focus on words that are representative of target domain to streamline
the adaptation process.

We propose a novel task-independent in-domain pre-training approach for
adapting PLMs that increases domain fit by focusing on keywords in the target
domain, where keywords are defined as “a sequence of one or more words that
offers a compact representation of a document’s content” (Rose et al., 2010).
By applying token masking only to in-domain keywords, the meaningful
information in the target domain is more directly captured by the PLM. This is
in contrast to the classic pre-training strategy that randomly masks tokens
(Devlin et al., 2019), which may overlook domain-meaningful information, or
the in-domain pre-training methods that selectively mask tokens deemed
important given the downstream task (Gu et al., 2020, inter alia), which
require incorporating information from the downstream task into the pre-
training stage. We empirically show that our method offers a better
transmission of high-quality information from the target domain into PLMs,
yielding better generalizability for the downstream tasks.

The key contributions of this paper are:

(1)

We propose the first task-agnostic selective masking technique for domain
adaptation of PLMs that relies solely on in-domain keywords. In particular, we
first extract contextually-relevant keywords from each available document in
the target domain using KeyBERT Grootendorst (2020) and keep the most
frequently occurring keywords to be masked during the adaptation phase.

(2)

We evaluate our proposed strategy by measuring the performance of fine-tuned
PLMs in six different settings. We leverage three different datasets for text
classification from multiple domains: IMDB movie reviews Maas et al. (2011),
Amazon pet product reviews from
Kaggle,222https://www.kaggle.com/datasets/kashnitsky/exploring-transfer-
learning-for-nlp and PUBHEALTH (Kotonya and Toni, 2020). Our experiments show
that the classifiers trained on top of two PLMs—in our case, Bidirectional
Encoder Representations from Transformers (BERT) Base and Large (Vaswani et
al., 2017; Devlin et al., 2019)—that are adapted based on our suggested
approach outperform all baselines, including the fine-tuned BERT with no in-
domain adaptation, and fine-tuned BERT adapted by random masking. Further, the
overhead of identifying in-domain keywords is reasonable, e.g., 7–15% of the
pre-training time (for two epochs of data) for BERT Large.

##  2 Related Work

Bidirectional Encoder Representations from Transformers (BERT) brought pre-
training to transformer networks Vaswani et al. (2017) through masked language
modeling (MLM) Devlin et al. (2019). They showed that a simple two-step
paradigm of generic pre-training followed by fine-tuning to the target domain
can significantly improve performance on a variety of tasks.

However, after showing that infusing an intermediate pre-training stage
(commonly known as in-domain pre-training) can help pre-trained Long Short-
Term Memory models learn domain-specific patterns better (Howard and Ruder,
2018; Peters et al., 2018), Gururangan et al. (2020) found that the same
advantage applies to PLMs as well. Since then, several efforts proposed
different domain-adaptive pre-training strategies.

Unsurprisingly, one of the most extensively utilized in-domain pre-training
methodologies has been to employ classic random masking to adapt PLMs into
several domains (Lee et al., 2020; Beltagy et al., 2019; Alsentzer et al.,
2019; Tavabi et al., 2022b, a; Araci, 2019). Following this, Zheng et al.
(2020) introduced the fully-explored MLM in which random masking is applied to
specific non-overlapping segments of the input sequence. The limitation of
random masking that we aim to address is that it may put unnecessary focus on
tokens that are not representative of the target domain.

In contrast, task-specific selective masking methods mask tokens that are
important to the downstream task. For each task, “importance“ is defined
differently: Gu et al. (2020) let an additional neural model learns important
tokens given the task at hand; Ziyadi et al. (2020) defined importance by
masking entities for the named entity recognition task, and Feng et al. (2018)
found important tokens by input reduction—maintaining model’s confidence in
the original prediction by reducing input—and they were left with a few
(potentially nonsensical) tokens that were treated as important to model.
Similarly, Li et al. (2020) designed a task-dependent objective for dialogue
adaptation, and Ke et al. (2019) proposed label-aware MLM for a sentiment
analysis task. In the same vein, token selection in certain domains, e.g.,
biomedical and clinical domains, was performed based on the entities relevant
to the domain Lin et al. (2021); Zhang et al. (2020b); Pergola et al. (2021).

Note that other MLM-based pre-training strategies focused on training a
language model from scratch (Zhang et al., 2020a; Joshi et al., 2020; Sun et
al., 2019, inter alia). However, since our work focuses on in-domain pre-
training, we skip this part for brevity.

In this study, we propose an information-based domain-adaptive pre-training
that, without being aware of the downstream task, selectively masks words that
are information-dense with respect to the target domain. As a result, PLMs
adapted using our mechanism outperform baselines adapted with random masking
or fine-tuned directly. In the following sections, we refer to our approach as
“keyword masking pre-training.“

##  3 Approach

###  3.1 Extracting In-domain Keywords

In order to extract keywords relevant to the domain of interest, we use
KeyBERT (Grootendorst, 2020). In a nutshell, KeyBERT uses BERT’s Devlin et al.
(2019) contextualized embeddings to find the $n$-grams–in our scenario,
unigrams–that concisely describe a given document. In particular, word
embeddings with the highest cosine similarity to the overall document-level
representation are identified as keywords that best represent the entire
document. We configure KeyBERT to extract up to 10 keywords from each input
document. Note that we do not pre-train or fine-tune BERT as the underlying
model for KeyBERT.

###  3.2 Removing Noisy Keywords

After extracting domain-specific keywords, we compute the frequency of each
specific word that has been recognized as a keyword in all in-domain
documents. Subsequently, we sort them in descending order of their frequency
and keep only the most frequent ones. This simple strategy allows us to remove
keywords that are likely to be noisy or irrelevant to the target domain.

Figure 1 summarizes the noisy keyword removal process for PUBHEALTH dataset
(see Appendix B for other domains). Note that the actual figure has a very
long tail on the right, indicating that the actual in-domain keywords (or
parts where information is condensed in the target domain) are frequently
repeated. The graph displays the frequency of terms along with the number of
times they are identified as keywords. In the PUBHEALTH dataset, for example,
more than 10,000 words were detected as keywords only once. Thus, we select
the cut-off point where the curve is intended to leap up, signaling that
keywords with repetition below the threshold were excluded from the list of
domain-relevant keywords.333The threshold is adjusted via three points: an
empirically chosen point from the graph, a point before, and a point after it.
Following keyword masking based on each of these three thresholds, we choose
the one that resulted in the highest F1 score on the validation split as the
final threshold. Namely, in the PUBHEALTH dataset, all words detected fewer
than eight times as a keyword were removed from the list of in-domain
keywords, and consequently, for performing keyword masking pre-training. The
provided examples on the graph in Figure 1 indicate a qualitative indication
that KeyBERT, coupled with our frequency-based heuristic, selects meaningful
domain-specific keywords. For example, our approach identifies relevant
keywords (e.g., health, coronavirus), while skipping other less relevant ones
(e.g., gym, gift).

Figure 1: The graph shows the frequency of the last 50 most frequent keywords
in the PUBHEALTH domain and the cut-off line for removing noisy keywords.

###  3.3 Keyword Masking Pre-training

We pair the list of retrieved candidate keywords with all target domain
documents to perform keyword masking pre-training. If any of the keywords from
the list appear in the input documents, the tokens corresponding to those
keywords get masked given the masking probability. In our pre-training
strategy, we use a constant learning rate scheduler with a high masking
probability rather than a linear one to force the majority of tokens
associated with keywords to be masked while continuously learning from
surrounding tokens. As our approach inherits from MLM (Devlin et al., 2019),
the tokens related to keywords are masked 80% of the time, replaced 10% of the
time with other tokens, and left unchanged 10% of the time. Note that during
pre-training masking only applies to the tokens that match the candidate
keywords. Therefore, there is no pre-training for unmasked tokens.

###  3.4 Fine-tuning and Baselines

We compare the performance of all fine-tuned PLMs adapted using our technique
with two other baselines: fine-tuned PLMs adapted using random masking, and
fine-tuned PLMs with no in-domain adaptation. For all these settings we employ
both BERT Base and BERT Large Devlin et al. (2019).

##  4 Experimental Setup

Data: In our experiments, we chose tasks and datasets with sufficient amounts
of unlabeled data for the domain adaptation stage in order to observe the
effects of keyword selection.444For example, we did not use the GLUE dataset
Wang et al. (2018) because the included texts are short. In particular, we
evaluate our method on three text classification datasets: PUBHEALTH Kotonya
and Toni (2020), which contains public health claims associated with veracity
labels, IMDB movie reviews dataset Maas et al. (2011), and Amazon pet product
reviews dataset (from a Kaggle
competition).555https://www.kaggle.com/datasets/kashnitsky/exploring-transfer-
learning-for-nlp

Based on the thresholds we studied for filtering out the noisy keywords (see
Section 3.2), we gathered 2,116, 7,274, and 6,881 domain-specific keywords
from the PUBHEALTH dataset, IMDB dataset, and Amazon dataset, respectively.

Settings:

We use KeyBERT Grootendorst (2020) to extract up to 10 unigram keywords per
input document utilizing contextualized word embeddings of BERT Base Devlin et
al. (2019), stratified by the Maximal Marginal Relevance (MMR) Carbonell and
Goldstein (1998) with a threshold of 0.8.

To perform keyword masking pre-training, we set the masking probability to
0.75 with a constant learning scheduler. The other hyperparameters are left at
their default values from the Hugging Face data collator for whole word
masking Wolf et al. (2020). For random masking pre-training, we set the
masking probability to 0.15, which is a standard value for continual MLM pre-
training, and left the remaining hyperparameters at the values provided by the
Hugging Face data collator for language modeling Wolf et al. (2020). Note that
the default learning rate scheduler is linear. Further, in all settings, pre-
training is limited to two epochs, and the batch size of 16 is adopted during
both the adaptation and fine-tuning stages.

With the learning rate set to 2e-5 and the weight decay set to 0.01 Devlin et
al. (2019), we fine-tune the whole network for all of our adapted models and
baselines for up to four epochs in all datasets, while keeping the other
hyperparameters at the default value of Hugging Face Wolf et al. (2020). The
models that obtained the highest F1 score in the validation partition are then
chosen and evaluated on the test split of the datasets.

_PUBHEALTH Dataset_
---
Adaptation Method | Accuracy (%) | F1 Score (%)
No Adaptation | 64.80 | 63.23
Random Masking | 65.77 | 64.94
Our Keyword Masking | *66.09 | *65.40
_IMDB Movie Reviews Dataset_
---
Adaptation Method | Accuracy (%) | F1 Score (%)
No Adaptation | 94.44 | 94.43
Random Masking | 94.96 | 94.95
Our Keyword Masking | *95.36 | *95.35
_Amazon Pet Product Reviews Dataset_
---
Adaptation Method | Accuracy (%) | F1 Score (%)
No Adaptation | 85.89 | 85.73
Random Masking | 86.33 | 86.31
Our Keyword Masking | *87.14 | *86.98
Table 1: A comparison between the performance of fine-tuning adapted PLMs
using our keyword masking and other baselines when BERT Base is used as the
PLM. The best results are shown bold and the obtained statistically
significant results compared to random masking are indicated by an asterisk
(*) (see Appendix E). _PUBHEALTH Dataset_
---
Adaptation Method | Accuracy (%) | F1 Score (%)
No Adaptation | 66.42 | 65.08
Random Masking | 63.90 | 64.74
Our Keyword Masking | *66.66 | 64.74
_IMDB Movie Reviews Dataset_
---
Adaptation Method | Accuracy (%) | F1 Score (%)
No Adaptation | 95.38 | 95.37
Random Masking | 95.50 | 95.49
Our Keyword Masking | 95.52 | 95.51
_Amazon Pet Product Reviews Dataset_
---
Adaptation Method | Accuracy (%) | F1 Score (%)
No Adaptation | 85.69 | 85.71
Random Masking | 86.84 | 86.72
Our Keyword Masking | *87.58 | *87.51
Table 2: A comparison between the performance of fine-tuning adapted PLMs
using our keyword masking and other baselines when BERT Large is used as the
PLM. The best results are shown bold and the obtained statistically
significant results compared to random masking are indicated by an asterisk
(*) (see Appendix E).

##  5 Results and Discussion

Table 1 and 2 report the performance of fine-tuned models that used multiple
domain-adaptive pre-training methods for each of our settings: three different
datasets and two distinct PLMs. Table 1 contains the results for BERT Base as
underlying PLM; Table 2 uses BERT Large.

In particular, each table contrasts the performance of two fine-tuned
baselines—one without adaptation/in-domain pre-training and one with random
masking in-domain pre-training—to a fine-tuned model adapted using our keyword
masking.

Both tables show that our approach outperforms all other baselines in all six
settings. The improvements are statistically significant in four out of six
settings (Appendix E). This highlights the importance of selecting
information-carrying keywords for masking during the in-domain pre-training.

The results reveal that our suggested in-domain pre-training technique
outperforms alternative settings with or without standard in-domain pre-
training on target domain unlabeled data. Although the benefits of continual
pre-training vary depending on the domain and the task at hand (Gururangan et
al., 2020), our adaptation strategy always has a greater impact on PLMs in
capturing domain-specific patterns compared to typical random masking when in-
domain adaptation has a positive impact on downstream tasks. This indicates
that our pre-training method indeed exposes the PLMs to relevant in-domain
representations.

Given the superior outcomes seen in our six different experiments, we can
argue that our selective masking strategy, which is task-agnostic as random
masking yet more effective, could potentially widely replace random masking in
the intermediate pre-training stage for a variety of NLP tasks. Other than
performance, our method is simple and has no “pathological behavior“ (Feng et
al., 2018) (see Appendix C). Additionally, our method takes 2 to 10 minutes of
computational overhead to extract keywords. This accounts for 7% to 39% of
pre-training time of only two epochs (Appendix A).

##  6 Conclusion

We proposed the first task-agnostic selective masking pre-training approach,
dubbed “keyword masking,“ to adapt PLMs to the target domains. For keyword
masking, we first extract in-domain keywords from the target domain using
KeyBERT (Grootendorst, 2020), and after excluding the noisy ones, we only mask
the selected keywords during adaptation.

We evaluated our methodology using six different settings. The results
revealed that when in-domain pre-training is conducted using our approach, all
fine-tuned PLMs outperform those with no adaptation or adapted using random
masking. Further, we observed that our pre-training approach was superior for
difficult tasks, i.e., datasets with many labels and more complexity. Lastly,
keyword masking pre-training can be widely substituted with random masking
during shift domain in NLP tasks since it is task-independent, as simple to
use as random masking, and more effective.

##  7 Limitations

Although all pre-training approaches require a sufficient amount of data,
given how we defined keywords, longer sequences suit our approach better than
short ones for studying the effects of keyword selection. Further, as shown in
this study, our findings strongly imply that the strategy we suggested for
adapting PLMs can effectively enhance their performance on text classification
as the downstream task. To determine whether these findings can translate to
other NLP applications, however, further experiments are required.

##  8 Ethics Statement

Although keyword extraction may amplify bias depending on the input documents
and the way it extracts keywords, KeyBERT (Grootendorst, 2020) has not been
reported to exhibit this behavior. Further work may be necessary to thoroughly
explore the potential of introducing undesired bias.

## References

  * Alsentzer et al. (2019) Emily Alsentzer, John R Murphy, Willie Boag, Wei-Hung Weng, Di Jin, Tristan Naumann, and Matthew McDermott. 2019.  Publicly available clinical bert embeddings.  _arXiv preprint arXiv:1904.03323_.
  * Araci (2019) Dogu Araci. 2019.  Finbert: Financial sentiment analysis with pre-trained language models.  _arXiv preprint arXiv:1908.10063_.
  * Beltagy et al. (2019) Iz Beltagy, Kyle Lo, and Arman Cohan. 2019.  Scibert: A pretrained language model for scientific text.  _arXiv preprint arXiv:1903.10676_.
  * Carbonell and Goldstein (1998) Jaime Carbonell and Jade Goldstein. 1998.  The use of mmr, diversity-based reranking for reordering documents and producing summaries.  In _Proceedings of the 21st Annual International ACM SIGIR Conference on Research and Development in Information Retrieval_ , SIGIR ’98, page 335–336, New York, NY, USA. Association for Computing Machinery.
  * Cohen (1960) Jacob Cohen. 1960.  A coefficient of agreement for nominal scales.  _Educational and Psychological Measurement_ , 20:37 – 46.
  * Devlin et al. (2019) Jacob Devlin, Ming-Wei Chang, Kenton Lee, and Kristina Toutanova. 2019.  BERT: Pre-training of deep bidirectional transformers for language understanding.  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 4171–4186, Minneapolis, Minnesota. Association for Computational Linguistics.
  * Efron (1979) B. Efron. 1979.  Bootstrap Methods: Another Look at the Jackknife.  _The Annals of Statistics_ , 7(1):1 – 26.
  * Efron (2003) Bradley Efron. 2003.  Second Thoughts on the Bootstrap.  _Statistical Science_ , 18(2):135 – 140.
  * Efron and Tibshirani (1993) Bradley Efron and Robert J. Tibshirani. 1993.  _An Introduction to the Bootstrap_.  Number 57 in Monographs on Statistics and Applied Probability. Chapman & Hall/CRC, Boca Raton, Florida, USA.
  * Feng et al. (2018) Shi Feng, Eric Wallace, Alvin Grissom II, Mohit Iyyer, Pedro Rodriguez, and Jordan Boyd-Graber. 2018.  Pathologies of neural models make interpretations difficult.  _arXiv preprint arXiv:1804.07781_.
  * Galton (1892) Francis Galton. 1892.  _Finger prints_.  57490-57492. Macmillan and Company.
  * Grootendorst (2020) Maarten Grootendorst. 2020.  Keybert: Minimal keyword extraction with bert.
  * Gu et al. (2020) Yuxian Gu, Zhengyan Zhang, Xiaozhi Wang, Zhiyuan Liu, and Maosong Sun. 2020.  Train no evil: Selective masking for task-guided pre-training.  In _Conference on Empirical Methods in Natural Language Processing_.
  * Gururangan et al. (2020) Suchin Gururangan, Ana Marasović, Swabha Swayamdipta, Kyle Lo, Iz Beltagy, Doug Downey, and Noah A. Smith. 2020.  Don’t stop pretraining: Adapt language models to domains and tasks.  _ArXiv_ , abs/2004.10964.
  * Howard and Ruder (2018) Jeremy Howard and Sebastian Ruder. 2018.  Universal language model fine-tuning for text classification.  _arXiv preprint arXiv:1801.06146_.
  * Joshi et al. (2020) Mandar Joshi, Danqi Chen, Yinhan Liu, Daniel S Weld, Luke Zettlemoyer, and Omer Levy. 2020.  Spanbert: Improving pre-training by representing and predicting spans.  _Transactions of the Association for Computational Linguistics_ , 8:64–77.
  * Ke et al. (2019) Pei Ke, Haozhe Ji, Siyang Liu, Xiaoyan Zhu, and Minlie Huang. 2019.  Sentilare: Sentiment-aware language representation learning with linguistic knowledge.  _arXiv preprint arXiv:1911.02493_.
  * Kotonya and Toni (2020) Neema Kotonya and Francesca Toni. 2020.  Explainable automated fact-checking: A survey.  In _Proceedings of the 28th International Conference on Computational Linguistics_ , pages 5430–5443, Barcelona, Spain (Online). International Committee on Computational Linguistics.
  * Lee et al. (2020) Jinhyuk Lee, Wonjin Yoon, Sungdong Kim, Donghyeon Kim, Sunkyu Kim, Chan Ho So, and Jaewoo Kang. 2020.  Biobert: a pre-trained biomedical language representation model for biomedical text mining.  _Bioinformatics_ , 36(4):1234–1240.
  * Li et al. (2020) Junlong Li, Zhuosheng Zhang, Hai Zhao, Xi Zhou, and Xiang Zhou. 2020.  Task-specific objectives of pre-trained language models for dialogue adaptation.  _ArXiv_ , abs/2009.04984.
  * Lin et al. (2021) Chen Lin, Timothy Miller, Dmitriy Dligach, Steven Bethard, and Guergana K. Savova. 2021.  Entitybert: Entity-centric masking strategy for model pretraining for the clinical domain.  In _Workshop on Biomedical Natural Language Processing_.
  * Maas et al. (2011) Andrew L. Maas, Raymond E. Daly, Peter T. Pham, Dan Huang, Andrew Y. Ng, and Christopher Potts. 2011.  Learning word vectors for sentiment analysis.  In _Proceedings of the 49th Annual Meeting of the Association for Computational Linguistics: Human Language Technologies_ , pages 142–150, Portland, Oregon, USA. Association for Computational Linguistics.
  * McHugh (2012) Mary L. McHugh. 2012.  Interrater reliability: the kappa statistic.  _Biochemia Medica_ , 22:276 – 282.
  * Pergola et al. (2021) Gabriele Pergola, Elena Kochkina, Lin Gui, Maria Liakata, and Yulan He. 2021.  Boosting low-resource biomedical qa via entity-aware masking strategies.  _arXiv preprint arXiv:2102.08366_.
  * Peters et al. (2018) Matthew E. Peters, Mark Neumann, Mohit Iyyer, Matt Gardner, Christopher Clark, Kenton Lee, and Luke Zettlemoyer. 2018.  Deep contextualized word representations.  In _Proceedings of the 2018 Conference of the North American Chapter of the Association for Computational Linguistics: Human Language Technologies, Volume 1 (Long Papers)_ , pages 2227–2237, New Orleans, Louisiana. Association for Computational Linguistics.
  * Rose et al. (2010) Stuart Rose, Dave Engel, Nick Cramer, and Wendy Cowley. 2010.  Automatic keyword extraction from individual documents.  _Text mining: applications and theory_ , 1(1-20):10–1002.
  * Sculley et al. (2015) David Sculley, Gary Holt, Daniel Golovin, Eugene Davydov, Todd Phillips, Dietmar Ebner, Vinay Chaudhary, Michael Young, Jean-Francois Crespo, and Dan Dennison. 2015.  Hidden technical debt in machine learning systems.  _Advances in neural information processing systems_ , 28.
  * Shin et al. (2020) Taylor Shin, Yasaman Razeghi, Robert L Logan IV, Eric Wallace, and Sameer Singh. 2020.  Autoprompt: Eliciting knowledge from language models with automatically generated prompts.  _arXiv preprint arXiv:2010.15980_.
  * Smeeton (1985) Nigel C. Smeeton. 1985.  Early history of the kappa statistic.  _Biometrics_ , 41(3):795–795.
  * Sun et al. (2019) Yu Sun, Shuohuan Wang, Yukun Li, Shikun Feng, Xuyi Chen, Han Zhang, Xin Tian, Danxiang Zhu, Hao Tian, and Hua Wu. 2019.  Ernie: Enhanced representation through knowledge integration.  _arXiv preprint arXiv:1904.09223_.
  * Tavabi et al. (2022a) Nazgol Tavabi, James Pruneski, Shahriar Golchin, Mallika Singh, Ryan Sanborn, Benton Heyworth, Amir Kimia, and Ata Kiapour. 2022a.  Building large-scale registries from unstructured clinical notes using a low-resource natural language processing pipeline.  _medRxiv_.
  * Tavabi et al. (2022b) Nazgol Tavabi, Marium Raza, Mallika Singh, Shahriar Golchin, Harsev Singh, Grant D Hogue, and Ata M Kiapour. 2022b.  A natural language processing pipeline to study disparities in cannabis use and documentation among children and young adults a survey of 21 years of electronic health records.  _medRxiv_.
  * Tunstall et al. (2022) L. Tunstall, L. von Werra, and T. Wolf. 2022.  _Natural Language Processing with Transformers_.  O’Reilly Media.
  * 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.  _Advances in neural information processing systems_ , 30.
  * Wang et al. (2018) Alex Wang, Amanpreet Singh, Julian Michael, Felix Hill, Omer Levy, and Samuel R Bowman. 2018.  Glue: A multi-task benchmark and analysis platform for natural language understanding.  _arXiv preprint arXiv:1804.07461_.
  * Wolf et al. (2020) Thomas Wolf, Lysandre Debut, Victor Sanh, Julien Chaumond, Clement Delangue, Anthony Moi, Pierric Cistac, Tim Rault, Remi Louf, Morgan Funtowicz, Joe Davison, Sam Shleifer, Patrick von Platen, Clara Ma, Yacine Jernite, Julien Plu, Canwen Xu, Teven Le Scao, Sylvain Gugger, Mariama Drame, Quentin Lhoest, and Alexander Rush. 2020.  Transformers: State-of-the-art natural language processing.  In _Proceedings of the 2020 Conference on Empirical Methods in Natural Language Processing: System Demonstrations_ , pages 38–45, Online. Association for Computational Linguistics.
  * Zhang et al. (2020a) Jingqing Zhang, Yao Zhao, Mohammad Saleh, and Peter Liu. 2020a.  Pegasus: Pre-training with extracted gap-sentences for abstractive summarization.  In _International Conference on Machine Learning_ , pages 11328–11339. PMLR.
  * Zhang et al. (2020b) Ningyu Zhang, Qianghuai Jia, Kangping Yin, Liang Dong, Feng Gao, and Nengwei Hua. 2020b.  Conceptualized representation learning for chinese biomedical text mining.  _arXiv preprint arXiv:2008.10813_.
  * Zheng et al. (2020) Mingzhi Zheng, Dinghan Shen, Yelong Shen, Weizhu Chen, and Lin Xiao. 2020.  Improving self-supervised pre-training via a fully-explored masked language model.  _arXiv preprint arXiv:2010.06040_.
  * Ziyadi et al. (2020) Morteza Ziyadi, Yuting Sun, Abhishek Goswami, Jade Huang, and Weizhu Chen. 2020\.  Example-based named entity recognition.  _arXiv preprint arXiv:2008.10570_.

##  Appendix A Cost of In-domain Keyword Extraction

For all of our settings, we leverage a single NVIDIA RTX A6000 GPU. Depending
on the size of the datasets and the length of input documents, running KeyBERT
(Grootendorst, 2020) for extracting in-domain keywords adds additional
computation that can be different between 2 and 10 minutes in our settings.
The overhead time for keyword extraction and the in-domain pre-training time
for each of the settings are compared in Table 3. As can be noticed, the time
ratio for keyword extraction to pre-training time ranges from 7% to 15% in
settings using BERT Large, and 19% to 39% for settings with BERT Base, which
is reasonable. Note that when pre-training is performed for more epochs, this
ratio noticeably decreases. The reported ratios are based on only two epochs
of in-domain pre-training in our settings.

Dataset Name |  BERT Base Pre-train Time |  BERT Large Pre-train Time |  Keyword Extraction Time |  Time Ratio to BERT Base (%) |  Time Ratio to BERT Large (%)
---|---|---|---|---|---
PUBHEALTH | 4.35 | 11.22 | 1.71 | 39 | 15
IMDB | 29.98 | 79.97 | 9.14 | 30 | 11
Amazon | 38.21 | 100.01 | 7.47 | 19 | 7
Table 3: The pre-training time for two epochs, and inference time for KeyBERT
in minutes.

##  Appendix B Removing Noisy Keywords (Graphs)

Similar to Figure 1, which illustrates the removal of noisy keywords for the
PUBHEALTH dataset, Figure 2 displays this procedure for the IMDB dataset and
the Amazon dataset.

Figure 2: The graphs show the frequency of the last 50 most frequent keywords
in IMDB and Amazon datasets along with the cut-off line for removing noisy
keywords. For keyword masking, keywords are selected from the subset to the
right of the cut-off line (due to space constraints, we do not show the actual
lengthy right tail of the charts). A few examples of words that were and were
not selected as in-domain keywords given this heuristic are shown on the
graphs as well.

##  Appendix C Pathological-free Behavior

It is possible that tokens to be selected for masking are not associated with
the domain according to human experts, but they nevertheless yield better
downstream classifiers. For instance, Feng et al. (2018) demonstrated that
even when the model is left with a small number of tokens after input
reduction, it can still be confident in its predictions even though the left
tokens are meaningless. A similar phenomenon was reported for prompting Shin
et al. (2020). To show our masking method’s non-pathological behavior, we
asked two human annotators to annotate the domain relevance of 50 randomly-
chosen words that were selected for masking by the respective method.666The
annotators were two of the authors. The annotations were independent, i.e., no
annotator saw the decisions of the other. The names of the methods used to
generate the 50 words to annotate were hidden during annotation. The
annotations were performed using a three-point Likert scale: irrelevant,
moderately relevant, and relevant.

Table 4 reports the results of this experiment as well as the Kappa inter-
annotator agreement score Galton (1892); Cohen (1960); Smeeton (1985); McHugh
(2012). We draw two observations from this table. First, the agreement between
the two annotators is high—substantial or near perfect—which indicates that
this task is well-defined. Second, the annotators agreed that the number of
moderately- or fully-relevant words is much higher in the keyword-based
strategy than in the random masking method. This result further highlights
that our masking strategy is indeed relying on the identification of domain-
relevant keywords to mask rather than picking up artifacts of the entanglement
present in neural architectures Sculley et al. (2015).

Dataset Name / Masking Method |  No. Irrelevant Words |  No. Moderately Related Words |  No. Related Words |  Kappa Value |  Level of Agreement
---|---|---|---|---|---
PUBHEALTH / Random Masking | 32 | 6 | 8 | 0.84 | Near Perfect
PUBHEALTH / Keyword Masking | 14 | 10 | 16 | 0.70 | Substantial
Amazon / Random Masking | 40 | 5 | 3 | 0.87 | Near Perfect
Amazon / Keyword Masking | 11 | 17 | 14 | 0.72 | Substantial
IMDB / Random Masking | 42 | 0 | 4 | 0.65 | Substantial
IMDB / Keyword Masking | 11 | 7 | 24 | 0.73 | Substantial
Table 4: The results of measuring inter-rater reliability using Cohen’s kappa
coefficient for 50 randomly selected words/tokens for masking during in-domain
pre-training.

##  Appendix D Implementation of Keyword Masking

To implement our keyword masking strategy, we develop a new data collator by
subclassing the Hugging Face data collator for whole word masking Wolf et al.
(2020). Our data collator masks only the tokens according to a certain list of
keywords given a probability of masking. Note that our data collator inherits
from MLM (Devlin et al., 2019), and no other words or tokens are masked during
pre-training except keywords provided by the list.

##  Appendix E Statistical Analysis

We analyze the statistical significance of the obtained improvements using a
bootstrap resampling technique with 1,000 samples in the resampling process
Efron (1979); Efron and Tibshirani (1993); Efron (2003). The hypothesis that
we investigate is if the results achieved by keyword masking are better than
the random masking pre-training strategy. We implement two variants of this
hypothesis: one compares F1 scores, and the other compares accuracies.

Dataset Name / PLM Name |  F1 Score $p$-value |  Accuracy $p$-value
---|---|---
PUBHEALTH / BERT Base | 0.015 | 0.018
PUBHEALTH / BERT Large | 0.505 | 0.010
IMDB / BERT Base | 0.046 | 0.050
IMDB / BERT Large | 0.468 | 0.454
Amazon / BERT Base | 0.000 | 0.000
Amazon / BERT Large | 0.002 | 0.002
Table 5: The computed $p$-values for F1 score and accuracy for each of our
settings using bootstrap resampling with 1,000 samples.

Table 5 lists the results of this analysis. Overall, the table exhibits that
in situations when PLM benefits well from in-domain pre-training, the
difference between keyword masking and random masking is statistically
significant for both F1 and accuracy scores with $p$-values $\leq 0.05$. The
differences are not statistically significant in the two scenarios: the IMDB
dataset with BERT Large and PUBHEALTH dataset with BERT Large (only for F1
score). This validates our findings since when BERT Large was employed, the
results from keyword masking and random masking for the IMDB dataset are quite
similar and close to the fine-tuned vanilla PLM. Similarly, in the PUBHEALTH
dataset with BERT Large, keyword masking and random masking tie in the F1
score, making the difference in the F1 score not statistically significant.
These results further confirm that the benefits of in-domain pre-training vary
depending on the domain and the task at hand Gururangan et al. (2020);
however, when in-domain pre-training has a positive impact on performance and
causes significant improvement compared to non-adapted setting, our approach
outperforms random masking and yields statistically significant gains.

##  Appendix F Detailed Description of Datasets

PUBHEALTH Dataset The PUBHEALTH dataset is divided into three sections: train,
test, and validation. Samples in each partition are public health claims with
one of four veracity labels including false, unproven, true, or mixture. The
labels were assigned by domain experts based on an explanation that they
provided for every claim, available in a separate column. These explanations
serve as in-domain unstructured data for our use. 9,832 samples in the train
split, 1,225 samples in the validation split, and 1,235 samples in the test
split form our dataset after a few unlabeled samples were removed.777This
dataset contains a small number of claims that did not fall under any of the
four aforementioned veracity labels.

IMDB Movie Reviews Dataset The two portions of the IMDB dataset are labeled
and unlabeled reviews, each having 50,000 reviews. The train, validation, and
test splits are generated by dividing the labeled portion by 80%, 10%, and
10%, respectively. That is, 40,000 reviews are allotted to the train split and
5,000 each to the validation and test splits. The unlabeled 50,000 reviews are
used for pre-training.

Amazon Pet Product Reviews Dataset There are six different labels for reviews
in the Amazon pet product dataset used in the Kaggle competition: dogs, fish
aquatic pets, cats, birds, bunny rabbit central, and small animals. The
dataset contains four splits: train, test, validation, and unlabeled. However,
since the test split does not include labels, we create our own test split by
randomly choosing a portion of the train split that is equal in size to the
validation split. As a result, in our setting the validation and test splits
each includes 17,353 samples; the train split contains 34,704 samples. In
addition, there are 100,000 reviews without labels in the dataset’s unlabeled
portion that serve as pre-training data.