Slavnov and Gaudin–Korepin Formulas for Models without ${\rm U}(1)$ Symmetry

Slavnov and Gaudin–Korepin Formulas for Models
without $\boldsymbol{{\rm U}(1)}$ Symmetry: the Twisted XXX Chain

Samuel BELLIARD † and Rodrigo A. PIMENTA ‡§

S. Belliard and R.A. Pimenta

† Laboratoire de Physique Théorique et Modélisation (CNRS UMR 8089),
† Université de Cergy-Pontoise, F-95302 Cergy-Pontoise, France
<EMAIL_ADDRESS>
‡ Departamento de Física, Universidade Federal de São Carlos,
‡ Caixa Postal 676, CEP 13565-905, São Carlos, Brasil<EMAIL_ADDRESS>
§ Physics Department, University of Miami, P.O. Box 248046, FL 33124, Coral
Gables, USA

Received September 02, 2015, in final form December 02, 2015; Published online
December 04, 2015

We consider the XXX spin-$\frac{1}{2}$ Heisenberg chain on the circle with an
arbitrary twist. We characterize its spectral problem using the modified
algebraic Bethe anstaz and study the scalar product between the Bethe vector
and its dual. We obtain modified Slavnov and Gaudin–Korepin formulas for the
model. Thus we provide a first example of such formulas for quantum integrable
models without ${\rm U}(1)$ symmetry characterized by an inhomogenous Baxter
T-Q equation.

algebraic Bethe ansatz; integrable spin chain; scalar product

82B23; 81R12

The study of quantum integrable models with ${\rm U}(1)$ symmetry by the Bethe
ansatz (BA) methods [4, 11, 34] provides exact solutions which found
applications in a wide range of domains such as: statistical physics,
condensed matter physics, high energy physics, mathematical physics, etc. One
of the major accomplishments of the method has been the obtaining of form
factors, for models related to ${\mathfrak{gl}}_{2}$ and ${\mathfrak{gl}}_{3}$
families of symmetry, in the compact form of determinants [8, 9, 21, 22, 24,
26, 31]. In particular, for models related to the ${\mathfrak{gl}}_{2}$
symmetry, the key results are the Slavnov [35] and the Gaudin–Korepin [19, 20,
25] formulas, which provide, respectively, the scalar product between an
eigenstate and an arbitrary state and the norm of the eigenstates.

In the case of models without ${\rm U}(1)$ symmetry, the usual BA techniques
in general fail to provide a complete description of the spectrum111In some
cases, some gauge transformation can allow to apply the ABA, see for example
the XYZ spin chain [4, 36]. For the XXX case, that we consider here, the ${\rm
GL}(2)$ symmetry allows one to restore the ${\rm U}(1)$ symmetry [14] and the
usual ABA applies, provided that the twist is a non-singular matrix. Also, in
the context of open XXZ spin chains, constraints on the parameters of the
model allow one to apply the usual techniques, see [28].. Thus alternative
methods have been developed, for instance, the separation of variables (SoV)
[15, 16, 23, 30, 33], the commuting transfer matrices method [3], the
functional method [18] or the $q$-Onsager approach [2]. Recently, key steps
have been accomplished for the Bethe ansatz solution of such models. On the
one hand, a new family of inhomogeneous Baxter T-Q equation to determine the
eigenvalues has been proposed by the off-diagonal Bethe ansatz (ODBA) [12,
37]. On the other hand, the construction of the off-shell Bethe vector has
been done in the context of a modified algebraic Bethe ansatz (MABA) approach
[1, 5, 6, 10]. Let us remember that previous developments in the BA technique,
in particular the obtainment of the eigenvectors of the XXX chain on the
segment with upper-triangular boundaries [7], brought important insights to
the MABA.

Here, we consider the question of the scalar product between the Bethe vectors
obtained from the modified algebraic Bethe ansatz. For simplicity, we consider
the case of the isotropic spin$-\frac{1}{2}$ Heisenberg chain on the circle
with an arbitrary twist. This is the simplest model which can be considered by
the MABA. In fact, its solution contains the main features entailed by the
method: the spectrum is characterized by an inhomogeneous Baxter T-Q relation
and the off-shell Bethe vector is generated by a modified creation operator.
In this context, we obtain a modified Slavnov formula for the scalar product
between an on-shell Bethe state and its off-shell dual. The formula (see (17))
is given in terms of a determinant depending on the Jacobian of the
inhomogenous eigenvalue expression. Moreover, it contains a new factor related
to a certain expansion of the Bethe vector. The square of the norm, i.e., the
modified Gaudin–Korepin formula (see (19)), is obtained by a limit. As
expected, the case with general integrable open boundary, which also breaks
the ${\rm U}(1)$ symmetry, has the same structure and will be presented in a
separated publication.

The isotropic spin$-\frac{1}{2}$ Heisenberg chain on the circle with an
arbitrary twist is given by the Hamiltonian

$\displaystyle
H=\,\sum_{k=1}^{N}\big{(}\sigma^{x}_{k}\otimes\sigma^{x}_{k+1}+\sigma^{y}_{k}\otimes\sigma^{y}_{k+1}+\sigma^{z}_{k}\otimes\sigma^{z}_{k+1}\big{)}$
(1)

subject to the following boundary conditions

$\displaystyle\gamma\sigma^{x}_{N+1}=\frac{\tilde{\kappa}^{2}+\kappa^{2}-\kappa_{+}^{2}-\kappa_{-}^{2}}{2}\sigma^{x}_{1}+i\frac{\kappa^{2}-\tilde{\kappa}^{2}-\kappa_{+}^{2}+\kappa_{-}^{2}}{2}\sigma^{y}_{1}+(\kappa\kappa_{-}-\tilde{\kappa}\kappa_{+})\sigma^{z}_{1},$
$\displaystyle\gamma\sigma^{y}_{N+1}=i\frac{\tilde{\kappa}^{2}-\kappa^{2}-\kappa_{+}^{2}+\kappa_{-}^{2}}{2}\sigma^{x}_{1}+\frac{\tilde{\kappa}^{2}+\kappa^{2}+\kappa_{+}^{2}+\kappa_{-}^{2}}{2}\sigma^{y}_{1}-i(\tilde{\kappa}\kappa_{+}+\kappa\kappa_{-})\sigma^{z}_{1},$
$\displaystyle\gamma\sigma^{z}_{N+1}=(\kappa\kappa_{+}-\tilde{\kappa}\kappa_{-})\sigma^{x}_{1}+i(\tilde{\kappa}\kappa_{-}+\kappa\kappa_{+})\sigma^{y}_{1}+(\tilde{\kappa}\kappa+\kappa_{+}\kappa_{-})\sigma^{z}_{1}.$

The twist parameters
$\\{\kappa,\tilde{\kappa},\kappa_{+},\kappa_{-}\\}\in\mathbb{C}^{4}$ are
generic and $\gamma=\tilde{\kappa}\kappa-\kappa_{+}\kappa_{-}$. The Pauli
matrices222$\sigma^{z}=\left(\begin{matrix}1&0\\\ 0&-1\end{matrix}\right)$,
$\sigma^{+}=\left(\begin{matrix}0&1\\\ 0&0\end{matrix}\right)$,
$\sigma^{-}=\left(\begin{matrix}0&0\\\ 1&0\end{matrix}\right)$,
$\sigma^{x}=\sigma^{+}+\sigma^{-}$, $\sigma^{y}=i(\sigma^{-}-\sigma^{+})$.
$\sigma^{\alpha}_{k}$ with $\alpha=x,y,z$ act non trivially on the $k$th space
of the quantum space ${\mathcal{H}}=\otimes_{k=1}^{N}V_{k}$ with
$V_{k}=\mathbb{C}^{2}$.

The Hamiltonian (1) is integrable and can be considered within the quantum
inverse scattering method [17]. Let us briefly recall this formalism. The key
object is the rational $R$-matrix

$\displaystyle R(u)=\frac{u}{c}+P,$

which acts on $\mathbb{C}^{2}\otimes\mathbb{C}^{2}$, with
$P=\sigma^{+}\otimes\sigma^{-}+\sigma^{-}\otimes\sigma^{+}+\frac{1}{2}(1+\sigma^{z}\otimes\sigma^{z})$
and $c\in\mathbb{C}^{*}$. From the $R$-matrix, we construct the monodromy
matrix

$\displaystyle T_{a}(u)=R_{a1}(u-\theta_{1})\cdots
R_{aN}(u-\theta_{N})=\left(\begin{matrix}t_{11}(u)&t_{12}(u)\\\
t_{21}(u)&t_{22}(u)\end{matrix}\right)_{a},$

which acts on $\mathbb{C}^{2}\otimes{\mathcal{H}}$ and with
$\\{\bar{\theta}\\}=\\{\theta_{1},\dots\theta_{N}\\}\in\mathbb{C}^{N}$ being
the inhomogeneity parameters. This monodromy matrix satisfies the RTT relation

$\displaystyle R_{ab}(u-v)T_{a}(u)T_{b}(v)=T_{b}(v)T_{a}(u)R_{ab}(u-v)$

that encodes commutation relations between the operators $\\{t_{ij}(u)\\}$,
see Appendix A. The transfer matrix, generating function of the conserved
quantities of the model, is given by

$\displaystyle
t(u)=\operatorname{Tr}_{a}\big{(}K_{a}T_{a}(u)\big{)}=\tilde{\kappa}t_{11}(u)+\kappa
t_{22}(u)+\kappa^{+}t_{21}(u)+\kappa^{-}t_{12}(u),$

with

$\displaystyle K=\left(\begin{matrix}\tilde{\kappa}&\kappa^{+}\\\
\kappa^{-}&\kappa\end{matrix}\right).$ (2)

The commutation relation between transfer matrices with different spectral
parameters

$\displaystyle[t(u),t(v)]=0,$

follows from the RTT relation and the ${\rm GL}(2)$ invariance of the
$R$-matrix

$\displaystyle[R_{ab}(u-v),K_{a}K_{b}]=0.$ (3)

The Hamiltonian (1) is given by

$\displaystyle H=2c\frac{d}{du}\big{(}\ln(t(u))\big{)}\big{|}_{u\to
0,\,\theta_{i}\to 0}-N,$

and thus its spectral problem is the same of the one of the transfer matrix.

The diagonalization of the transfer matrix can be obtained by means of the
MABA [1, 5, 6, 10] and leads to an inhomogeneous Baxter T-Q equation. In order
to do that, we introduce the following transformation of the twist matrix (2)

$\displaystyle K=LDL$

with

$\displaystyle
L=\mu^{\frac{1}{2}}\left(\begin{matrix}1&\frac{\rho}{\kappa^{-}}\\\
\frac{\rho}{\kappa^{+}}&1\end{matrix}\right),\qquad
D=\left(\begin{matrix}\tilde{\kappa}-\rho&0\\\
0&\kappa-\rho\end{matrix}\right),\qquad\mu=\frac{\tilde{\kappa}+\kappa-\rho}{\tilde{\kappa}+\kappa-2\rho}$

and

$\displaystyle\rho^{2}-(\tilde{\kappa}+\kappa)\rho+\kappa^{+}\kappa^{-}=0.$

By means of this transformation, we can obtain a modified monodromy matrix
$\overline{T}(u)=LT(u)L$ with entries given by modified operators
$\\{\nu_{ij}(u)\\}$. They are expressed in terms of the initial
$\\{t_{ij}(u)\\}$ operators by

$\displaystyle\nu_{11}(u)=\mu\left(t_{11}(u)+\frac{\rho}{\kappa^{+}}t_{12}(u)+\frac{\rho}{\kappa^{-}}t_{21}(u)+\frac{\rho^{2}}{\kappa^{+}\kappa^{-}}t_{22}(u)\right),$
$\displaystyle\nu_{12}(u)=\mu\left(t_{12}(u)+\frac{\rho}{\kappa^{-}}(t_{11}(u)+t_{22}(u))+\left(\frac{\rho}{\kappa^{-}}\right)^{2}t_{21}(u)\right),$
$\displaystyle\nu_{21}(u)=\mu\left(t_{21}(u)+\frac{\rho}{\kappa^{+}}(t_{11}(u)+t_{22}(u))+\left(\frac{\rho}{\kappa^{+}}\right)^{2}t_{12}(u)\right),$
$\displaystyle\nu_{22}(u)=\mu\left(t_{22}(u)+\frac{\rho}{\kappa^{+}}t_{12}(u)+\frac{\rho}{\kappa^{-}}t_{21}(u)+\frac{\rho^{2}}{\kappa^{+}\kappa^{-}}t_{11}(u)\right).$

It follows that the transfer matrix has a modified diagonal form given by

$\displaystyle
t(u)=\operatorname{Tr}\big{(}D\bar{T}(u)\big{)}=(\tilde{\kappa}-\rho)\nu_{11}(u)+(\kappa-\rho)\nu_{22}(u).$

To construct the Bethe vector we use the usual highest weight representation
and the highest weight vector

$\displaystyle|0\rangle=\left(\begin{matrix}1\\\
0\end{matrix}\right)^{\otimes^{N}}.$

The actions of the operators $\\{t_{ij}(u)\\}$ on it are given by

$\displaystyle t_{ii}(u)|0\rangle=\lambda_{i}(u)|0\rangle,\qquad
t_{21}(u)|0\rangle=0,$

with

$\displaystyle\lambda_{1}(u)=\prod_{i=1}^{N}\frac{u-\theta_{i}+c}{c},\qquad\lambda_{2}(u)=\prod_{i=1}^{N}\frac{u-\theta_{i}}{c}.$
(4)

For the dual Bethe vector we use the dual highest weight vector

$\displaystyle\langle 0|=(1,\,0)^{\otimes^{N}},$

with the actions

$\displaystyle\langle 0|t_{ii}(u)=\lambda_{i}(u)\langle 0|,\qquad\langle
0|t_{12}(u)=0,\qquad\langle 0|0\rangle=1.$

For the operators $\\{\nu_{ij}(u)\\}$, we have a modified action on the
highest weight vector. We can show that

$\displaystyle\nu_{11}(u)|0\rangle=\lambda_{1}(u)|0\rangle+\frac{\rho}{\kappa^{+}}\nu_{12}(u)|0\rangle,$
(5)
$\displaystyle\nu_{22}(u)|0\rangle=\lambda_{2}(u)|0\rangle+\frac{\rho}{\kappa^{+}}\nu_{12}(u)|0\rangle,$
(6)
$\displaystyle\nu_{21}(u)|0\rangle=\frac{\rho}{\kappa^{+}}\big{(}\lambda_{1}(u)+\lambda_{2}(u)\big{)}|0\rangle+\left(\frac{\rho}{\kappa^{+}}\right)^{2}\nu_{12}(u)|0\rangle.$
(7)

We have thus all the ingredients to implement the MABA. We will use the
notation $\bar{u}$ with $\\#\bar{u}=M$ for the set of $M$ variables
$\\{u_{1},u_{2},\dots,u_{M}\\}$. If the element $u_{i}$ is removed, we note
$\bar{u}_{i}=\\{u_{1},u_{2},\dots,u_{i-1},u_{i+1},\dots,u_{M}\\}$. If we also
remove the element $u_{j}$, we note $\bar{u}_{ij}=\bar{u}/\\{u_{i},u_{j}\\}$.
For products of functions (see (8) bellow) or of operators
$\left\\{\nu_{ij}(u)\right\\}$, we use the convention

$\displaystyle g(u,\bar{u})=\prod_{i=1}^{M}g(u,u_{i}),\qquad
g(\bar{v},\bar{u})=\prod_{i=1}^{M}\prod_{j=1}^{M}g(v_{j},u_{i}),$
$\displaystyle g(u_{i},\bar{u}_{i})=\prod_{j=1,j\neq
i}^{M}g(u_{i},u_{j}),\qquad\nu_{ij}(\bar{u})=\prod_{k=1}^{M}\nu_{ij}(u_{k}).$

The functions

$\displaystyle g(u,v)=\frac{c}{u-v},\qquad f(u,v)=1+g(u,v)=\frac{u-v+c}{u-v}$
(8)

will be widely used.

Let us consider the vector

$\displaystyle\mathbb{B}^{M}(\bar{u})=\nu_{12}(u_{1})\cdots\nu_{12}(u_{M})|0\rangle=\nu_{12}(\bar{u})|0\rangle,$
(9)

and act with the transfer matrix on it. In order to perform this calculation
we need to find the action of the operators $\\{\nu_{ij}(u)\\}$ on the vector
(9). Using the ${\rm GL}(2)$ invariance (3), it is easy to see that the new
operators satisfy the same commutation relations (23), (24), (25) of the
operators $\\{t_{ij}(u)\\}$, see Appendix A. Thus, using these commutation
relations and the action (5), (6), (7), we can show that

$\displaystyle\nu_{12}(u)\mathbb{B}^{M}(\bar{u})=\mathbb{B}^{M+1}(u,\bar{u}),$
$\displaystyle\nu_{11}(u)\mathbb{B}^{M}(\bar{u})=\frac{\rho}{\kappa^{+}}\mathbb{B}^{M+1}(u,\bar{u})+\lambda_{1}(u)f(\bar{u},u)\mathbb{B}^{M}(\bar{u})$
$\displaystyle\hphantom{\nu_{11}(u)\mathbb{B}^{M}(\bar{u})=}{}+\sum_{i=1}^{M}g(u,u_{i})\lambda_{1}(u_{i})f(\bar{u}_{i},u_{i})\mathbb{B}^{M}(u,\bar{u}_{i}),$
$\displaystyle\nu_{22}(u)\mathbb{B}^{M}(\bar{u})=\frac{\rho}{\kappa^{+}}\mathbb{B}^{M+1}(u,\bar{u})+\lambda_{2}(u)f(u,\bar{u})\mathbb{B}^{M}(\bar{u})$
$\displaystyle\hphantom{\nu_{22}(u)\mathbb{B}^{M}(\bar{u})=}{}+\sum_{i=1}^{M}g(u_{i},u)\lambda_{2}(u_{i})f(u_{i},\bar{u}_{i})\mathbb{B}^{M}(u,\bar{u}_{i}),$
$\displaystyle\nu_{21}(u)\mathbb{B}^{M}(\bar{u})=\left(\frac{\rho}{\kappa^{+}}\right)^{2}\mathbb{B}^{M+1}(u,\bar{u})$
(10)
$\displaystyle\hphantom{\nu_{21}(u)\mathbb{B}^{M}(\bar{u})=}{}+\frac{\rho}{\kappa^{+}}\left(\Lambda^{M}_{d}(u,\bar{u}|1,1)\mathbb{B}^{M}(\bar{u})+\sum_{i=1}^{M}g(u_{i},u)E^{M}_{d}(u_{i},\bar{u}_{i}|1,1)\mathbb{B}^{M}(u,\bar{u}_{i})\right)$
$\displaystyle\hphantom{\nu_{21}(u)\mathbb{B}^{M}(\bar{u})=}{}+\left(\sum_{i=1}^{M}F(u,u_{i},\bar{u}_{i})\mathbb{B}^{M-1}(\bar{u}_{i})+\sum_{1\leq
i<j\leq
M}G(u,u_{i},u_{j},\bar{u}_{ij})\mathbb{B}^{M-1}(u,\bar{u}_{ij})\right),$

with

$\displaystyle\Lambda^{M}_{d}(u,\bar{u}|x,y)=xf(\bar{u},u)\lambda_{1}(u)+yf(u,\bar{u})\lambda_{2}(u),$
$\displaystyle
E^{M}_{d}(u_{i},\bar{u}_{i}|x,y)=-xf(\bar{u}_{i},u_{i})\lambda_{1}(u_{i})+yf(u_{i},\bar{u}_{i})\lambda_{2}(u_{i}),$
$\displaystyle
F(u,u_{i},\bar{u}_{i})=g(u,u_{i})\lambda_{1}(u)\lambda_{2}(u_{i})f(u,\bar{u}_{i})f(\bar{u}_{i},u_{i})+g(u_{i},u)\lambda_{1}(u_{i})\lambda_{2}(u)f(u_{i},\bar{u}_{i})f(\bar{u}_{i},u),$
$\displaystyle
G(u,u_{i},u_{j},\bar{u}_{ij})=g(u,u_{i})g(u_{j},u)\lambda_{1}(u_{i})\lambda_{2}(u_{j})f(u_{i},u_{j})f(u_{i},\bar{u}_{ij})f(\bar{u}_{ij},u_{j})$
$\displaystyle\hphantom{G(u,u_{i},u_{j},\bar{u}_{ij})=}{}+g(u,u_{j})g(u_{i},u)\lambda_{1}(u_{j})\lambda_{2}(u_{i})f(u_{j},u_{i})f(u_{j},\bar{u}_{ij})f(\bar{u}_{ij},u_{i}),$

and where we have used the functional identities

$\displaystyle
f(\bar{u},u)+\sum_{i=1}^{M}g(u,u_{i})f(\bar{u}_{i},u_{i})=1,\qquad
f(u,\bar{u})+\sum_{i=1}^{M}g(u_{i},u)f(u_{i},\bar{u}_{i})=1.$

It follows that the action of the transfer matrix is given by

$\displaystyle
t(u)\mathbb{B}^{M}(\bar{u})=\frac{\kappa^{-}}{\mu}\mathbb{B}^{M+1}(u,\bar{u})+\Lambda^{M}_{d}(u,\bar{u}|\tilde{\kappa}-\rho,\kappa-\rho)\mathbb{B}^{M}(\bar{u})$
$\displaystyle\hphantom{t(u)\mathbb{B}^{M}(\bar{u})=}{}+\sum_{i=1}^{M}g(u_{i},u)E^{M}_{d}(u_{i},\bar{u}_{i}|\tilde{\kappa}-\rho,\kappa-\rho)\mathbb{B}^{M}(u,\bar{u}_{i}),$
(11)

where we have used the relation
$\frac{\rho}{\kappa^{+}}(\tilde{\kappa}+\kappa-2\rho)=\frac{\kappa^{-}}{\mu}$.
The new term $\mathbb{B}^{M+1}(u,\bar{u})$ is characteristic for models which
break the ${\rm U}(1)$ symmetry. From the formula (11) we can obtain, by limit
the upper triangular case $\kappa^{-}=0$ and the diagonal case
$\kappa^{-}=\kappa^{+}=0$. The Bethe ansatz solution is then obtained by
requiring that
$E^{M}_{d}(u_{i},\bar{u}_{i}|\tilde{\kappa}-\rho,\kappa-\rho)=0$ for
$M=0,\dots,N$. Only the part with positive total spin contributes to the
solution. Here to complete the MABA we must thus find the action of the
operator $\nu_{12}(u)$ on the Bethe vector (9) when $M=N$. In this case the
$2^{N}$ independent diagonal Bethe vectors with $M\leq N$ and partitions of
$N$ variables are contained in the vector, see Appendix B for the expression
of the Bethe vector in terms of the initial operator $t_{12}(u)$. This allows
us to find an off-shell action with a wanted/unwanted form

$\displaystyle\frac{\kappa^{-}}{\mu}\mathbb{B}^{N+1}(u,\bar{u})=\Lambda_{g}^{N}(u,\bar{u})\mathbb{B}^{N}(\bar{u})+\sum_{i=1}^{N}g(u_{i},u)E_{g}^{N}(u_{i},\bar{u}_{i})\mathbb{B}^{N}(u,\bar{u}_{i}),$
(12)
$\displaystyle\Lambda_{g}^{N}(u,\bar{u})=2\rho\lambda_{1}(u)\lambda_{2}(u)g(u,\bar{u}),\qquad
E_{g}^{N}(u_{i},\bar{u}_{i})=2\rho\lambda_{1}(u_{i})\lambda_{2}(u_{i})g(u_{i},\bar{u}_{i}).$

This action can be proved following the method given in [13] and will be
considered elsewhere. Finally, from (11) and (12), we obtain the off-shell
equation with generic $\bar{u}$ with $\\#\bar{u}=N$,

$\displaystyle
t(u)\,\mathbb{B}^{N}(\bar{u})=\Lambda^{N}(u,\bar{u})\mathbb{B}^{N}(\bar{u})+\sum_{i=1}^{N}g(u_{i},u)E^{N}(u_{i},\bar{u}_{i})\mathbb{B}^{N}(u,\bar{u}_{i}),$

with an inhomogeneous eigenvalue

$\displaystyle\Lambda^{N}(u,\bar{u})=(\tilde{\kappa}-\rho)\lambda_{1}(u)f(\bar{u},u)+(\kappa-\rho)\lambda_{2}(u)f(u,\bar{u})+2\rho\lambda_{1}(u)\lambda_{2}(u)g(u,\bar{u}),$
(13)

and an inhomogeneous Bethe equation

$\displaystyle
E^{N}(u_{i},\bar{u}_{i})=-(\tilde{\kappa}-\rho)\lambda_{1}(u_{i})f(\bar{u}_{i},u_{i})+(\kappa-\rho)\lambda_{2}(u_{i})f(u_{i},\bar{u}_{i})$
$\displaystyle\hphantom{E^{N}(u_{i},\bar{u}_{i})=}{}+2\rho\lambda_{1}(u_{i})\lambda_{2}(u_{i})g(u_{i},\bar{u}_{i}),$
(14)

for $i=1,\dots,N$.

One can proceed in a similar way for the dual Bethe vector

$\displaystyle\mathbb{C}^{N}(\bar{u})=\langle 0|\nu_{21}(\bar{u}),$ (15)

and, in particular, obtain

$\displaystyle\mathbb{C}^{N}(\bar{u})t(u)=\Lambda^{N}(u,\bar{u})\mathbb{C}^{N}(\bar{u})+\sum_{i=1}^{N}g(u_{i},u)E^{N}(u_{i},\bar{u}_{i})\mathbb{C}^{N}(u,\bar{u}_{i}),$

with the same eigenvalue and Bethe equations of the Bethe vector (9) with
$M=N$.

When the Bethe equations are satisfied, i.e., $E^{N}(u_{i},\bar{u}_{i})=0$ for
$i=1,\dots,N$, and we consider non-singular solutions of the Bethe equations
[29], the on-shell Bethe vectors are eigenstates of the transfer matrix

$\displaystyle
t(u)\mathbb{B}^{N}(\bar{u})=\Lambda^{N}(u,\bar{u})\mathbb{B}^{N}(\bar{u}),\qquad\mathbb{C}^{M}(\bar{u})t(u)=\Lambda^{N}(u,\bar{u})\mathbb{C}^{N}(\bar{u}).$

Let us remark that the completeness of the solution given by (13), (14) has
been numerically checked for chains with small size. It should be interesting
to prove it along the lines of [27].

We are now in position to consider the scalar products for the Bethe vector
(9) and (15), namely,

$\displaystyle
S^{N}(\bar{u}|\bar{v})=\mathbb{C}^{N}(\bar{u})\mathbb{B}^{N}(\bar{v}).$ (16)

From the construction given hereafter, we obtain the modified Slavnov formula
of the twisted XXX spin chain characterized by the inhomogeneous Baxter T-Q
equation (13). When

$\displaystyle E^{N}(u_{i},\bar{u}_{i})=0,$

for $i=1,\dots,N$, the scalar product (16) has a compact form given by

$\displaystyle\hat{S}^{N}(\bar{u},\bar{v})=c^{N}\left(\frac{\mu^{2}}{\tilde{\kappa}+\kappa-\rho}\right)^{N}W_{0}^{N}(\bar{u})\frac{\operatorname{Det}_{N}\big{(}\frac{\partial}{\partial
u_{i}}\Lambda^{N}(v_{j},\bar{u})\big{)}}{\operatorname{Det}_{N}\big{(}g(v_{i},u_{j})\big{)}}$
(17)

with

$\displaystyle
W_{0}^{N}(\bar{u})=\left(\frac{\kappa^{-}}{\mu\rho}\right)^{N}\langle
0|\nu_{12}(\bar{u})|0\rangle,$

which is given explicitly in Appendix B. When

$\displaystyle E^{N}(v_{i},\bar{v}_{i})=0,$

for $i=1,\dots,N$, the scalar product (16) is given by

$\displaystyle\tilde{S}^{N}(\bar{u},\bar{v})=c^{N}\left(\frac{\mu^{2}}{\tilde{\kappa}+\kappa-\rho}\right)^{N}W_{0}^{N}(\bar{v})\frac{\operatorname{Det}_{N}\big{(}\frac{\partial}{\partial
v_{i}}\Lambda^{N}(u_{j},\bar{v})\big{)}}{\operatorname{Det}_{N}\big{(}g(u_{i},v_{j})\big{)}}.$
(18)

We can obtain the square of the norm by imposing the limit $\bar{u}=\bar{v}$.
Using the well-known formula for the Cauchy determinant

$\displaystyle\operatorname{Det}_{N}\big{(}g(v_{i},u_{j})\big{)}=\frac{g(\bar{v},\bar{u})}{\prod\limits_{i<j}g(u_{i},u_{j})g(v_{j},v_{i})},$

we obtain

$\displaystyle{\mathscr{N}}^{N}(\bar{u})=\left(\frac{\mu^{2}}{\tilde{\kappa}+\kappa-\rho}\right)^{N}W_{0}^{N}(\bar{u})\left(\prod_{i<j}^{N}g(u_{i},u_{j})g(u_{j},u_{i})\right)\operatorname{Det}_{N}(G_{ij}),$
(19)

where the matrix elements $G_{ij}$, for $i,j=1,\dots,N$, are given by

$\displaystyle G_{ii}=2\rho
c(\lambda_{2}(u_{i})\partial_{u_{i}}\lambda_{1}(u_{i})+\lambda_{1}(u_{i})\partial_{u_{i}}\lambda_{2}(u_{i}))$
$\displaystyle\hphantom{G_{ii}=}{}+(-1)^{N}(\tilde{\kappa}-\rho)\left(c\,h(\bar{u}_{i},u_{i})\partial_{u_{i}}\lambda_{1}(u_{i})-\lambda_{1}(u_{i})\sum_{j=1,j\neq
i}^{N}h(\bar{u}_{ij},u_{i})\right)$
$\displaystyle\hphantom{G_{ii}=}{}+(\kappa-\rho)\left(ch(u_{i},\bar{u}_{i})\partial_{u_{i}}\lambda_{2}(u_{i})+\lambda_{2}(u_{i})\sum_{j=1,j\neq
i}^{N}h(u_{i},\bar{u}_{ij})\right)$ $\displaystyle
G_{ij}=(-1)^{N}(\tilde{\kappa}-\rho)\lambda_{1}(u_{j})h(\bar{u}_{ij},u_{j})-(\kappa-\rho)\lambda_{2}(u_{j})h(u_{j},\bar{u}_{ij})\qquad\text{for}\quad
i\neq j,$

with

$\displaystyle h(u,v)=\frac{f(u,v)}{g(u,v)}=\frac{u-v+c}{c},$

and where we have applied the l’Hospital’s rule to find

$\displaystyle G_{ij}=\lim_{v_{j}\to u_{j}}c\frac{\frac{\partial}{\partial
u_{i}}\Lambda^{N}(v_{j},\bar{u})}{g(v_{j},\bar{u})}.$ (20)

From the conjectured modified Slavnov and Gaudin–Korepin formulas we remark
that the ratio, needed for the calculation of the correlations functions, is
independent of the $W_{0}^{N}(\bar{u})$ and of the constant
$\big{(}\frac{\mu^{2}}{\tilde{\kappa}+\kappa-\rho}\big{)}^{N}$. Thus the
relevant part of the formula for the correlations functions is given by the
Jacobian of the inhomogeneous eigenvalue (13) and it limits (20).

To obtain this conjecture we have proceeded in the following way. We start
from the following hypothesis:

  * •

We can impose the Bethe equations (14), $E^{N}(v_{i},\bar{v}_{i})=0$ for
$i=1,\dots,N$, by linearizing the quadratic terms
$\lambda_{1}(v_{i})\lambda_{2}(v_{i})$ in terms of $\lambda_{1}(v_{i})$ and
$\lambda_{2}(v_{i})$.

  * •

The usual Slavnov formula must be recovered in the ${\rm U}(1)$ symmetric
limit and thus the modified Slavnov formula contain the determinant of the
Jacobian of the inhomogenous eigenvalue (13).

For the ${\rm U}(1)$ symmetric case ($\kappa^{+}=\kappa^{-}=0$), the scalar
product between an on-shell Bethe vector and an off-shell Bethe vector is
given by the Slavnov formula [35]. If the $\bar{v}$ are on-shell, i.e.,
$E^{M}_{d}(v_{i},\bar{v}_{i}|\tilde{\kappa},\kappa)=0$ for $i=1,\dots,M$, we
have

$\displaystyle\tilde{S}_{d}^{M}(\bar{u},\bar{v})=\left(\frac{c}{\tilde{\kappa}}\right)^{M}\lambda_{2}(\bar{v})\frac{{\operatorname{Det}_{M}\big{(}\frac{\partial}{\partial
v_{i}}\Lambda_{d}^{M}(u_{j},\bar{v},\tilde{\kappa},\kappa)\big{)}}}{\operatorname{Det}_{M}(g(u_{i},v_{j}))}.$

From the first term of action the (10), for $M=N$ and using (12), we can
extract the leading term of the scalar product with $N$ quadratic terms:
$\lambda_{1}(u_{i})\lambda_{2}(u_{i})$ with $i=1,\dots,N$. If we impose the
Bethe equations (14), $E^{N}(v_{i},\bar{v}_{i})=0$ for $i=1,\dots,N$, by
linearizing the quadratic terms $\lambda_{1}(v_{i})\lambda_{2}(v_{i})$ in
terms of $\lambda_{1}(v_{i})$ and $\lambda_{2}(v_{i})$, this leading term is
invariant and is the only contribution at the top order $3N$ in the
$\lambda_{i}$. This term is proportional to

$\displaystyle
W_{0}^{N}(\bar{v})\prod_{i=1}^{N}\Lambda^{N}_{g}(u_{i},\bar{v}).$

It allows us to fix the functional to complete the determinant part of the
Slavnov formula (18). Indeed, the functional is the same, up to a constant,
when we consider the leading coefficient in the $\lambda_{i}$ and use the
identity

$\displaystyle\operatorname{Det}_{N}\left(\frac{\partial}{\partial
v_{i}}\Lambda_{g}^{N}(u_{j},\bar{v})\right)=\frac{1}{c^{N}}\operatorname{Det}_{N}\big{(}g(u_{j},v_{i})\big{)}\prod_{i=1}^{N}\Lambda_{g}^{N}(u_{i},\bar{v}).$

In the ${\rm U}(1)$ symmetry limit, $\rho=0$, we have
$W_{0}^{N}(\bar{v})=\lambda_{2}(\bar{v})(\frac{\kappa}{\tilde{\kappa}}+1)^{N}$
and restore the usual Slavnov formula up to a constant. The formula (17) can
de derived, also up to a constant, in a similar way starting from the action
of the operator $\nu_{12}(u)$ on the dual Bethe vector (15). This way to find
the modified Slavnov formula does not fix the constant.

Let us fix this constant by considering the simplest case $N=1$ and discuss
the difficulties to find the good parametrization for the off-shell scalar
product, which allows us to find the modified Slavnov formula by linearizing
the quadratic term of the Bethe equation, for general $N$. In the case $N=1$,
the good one is

$\displaystyle
S^{1}(u,v)=\mu\left(S_{d}^{1}(u,v)+\frac{\mu}{\tilde{\kappa}+\kappa-\rho}\big{(}\Lambda^{1}_{g}(u,v)W_{0}^{1}(v)+\Lambda^{1}_{g}(v,u)W_{0}^{1}(u)\big{)}\right)$
(21)

with

$\displaystyle W_{0}^{1}(u)=\lambda_{1}(u)+\lambda_{2}(u),\qquad
S_{d}^{1}(u,v)=g(u,v)(\lambda_{1}(v)\lambda_{2}(u)-\lambda_{1}(u)\lambda_{2}(v)).$

To obtain the Slavnov formula we use the Bethe equation to linearize the
quadratic term that depends on $v$

$\displaystyle\lambda_{1}(v)\lambda_{2}(v)=\frac{1}{2\rho}\big{(}(\tilde{\kappa}-\rho)\lambda_{1}(v)-(\kappa-\rho)\lambda_{2}(v)\big{)}$
(22)

and thus we arrive to our result

$\displaystyle
S_{E(v)=0}^{1}(u,v)=\frac{\mu^{2}}{\tilde{\kappa}+\kappa-\rho}g(v,u)W_{0}^{1}(v)\big{(}(\tilde{\kappa}-\rho)\lambda_{1}(u)-(\kappa-\rho)\lambda_{2}(u)-2\rho\lambda_{1}(u)\lambda_{2}(u)\big{)}.$

If we also linearize the second quadratic term, we find zero for $u\neq v$,
which shows the orthogonality of the Bethe vectors and the modified
Gaudin–Korepin formula (19) for $u=v$.

The parametrization (21) is not unique; for instance, from the projection on
the $\\{t_{ij}(u)\\}$ operator, see Appendix B, we find

$\displaystyle
S^{1}(u,v)=\mu^{2}S_{d}^{1}(u,v)+\mu^{2}\frac{\rho^{2}}{\kappa^{+}\kappa^{-}}\big{(}\lambda_{1}(u)+\lambda_{2}(u)\big{)}\big{(}\lambda_{1}(v)+\lambda_{1}(v)\big{)}$

that is equivalent to (21) when we specify the explicit form of the
$\lambda_{i}(u)$ (4). In this case the quadratic term
$\lambda_{1}(v)\lambda_{2}(v)$ does not appear and the formula (22) could not
be used directly. Moreover, the recursion relation (10) gives also another
parametrization in order 3 in the $\lambda_{i}(u)$ like (21) but with terms of
the form $\lambda_{i}(u)\lambda_{i}(v)$. In this case the linearization of the
quadratic term (22) could be used but does not lead directly to the modified
Slavnov formula. A systematic way to fix the good parametrization of the off-
shell scalar product (i.e., a parametrization that reduces to the modified
Slavnov formula, when we linearize the quadratic terms
$\lambda_{1}(v_{i})\lambda_{2}(v_{i})$ from the inhomogenous Bethe equations)
remains an open problem for the moment and will be discussed elsewhere. The
conjecture can be tested exactly for $N=1$, which allows us to fix the
constant, and then the cases $N=2$ and $N=3$ have been checked numerically to
support the conjecture.

Finally, let us point out another problem which is worth to be considered: the
construction of an explicit algebraic link between the MABA solution and the
solution obtained by means of the usual ABA, see, e.g., [32]. One first step
to address this problem, at least at the spectral level, is to equate the
eigenvalue expression found from the MABA (13) with the ones obtained from the
usual ABA. The simplest eigenvalue from the usual ABA is given by

$\displaystyle\alpha\lambda_{1}(u)+(\kappa+\tilde{\kappa}-\alpha)\lambda_{2}(u)$

with
$\alpha=\frac{1}{2}(\kappa+\tilde{\kappa}+\sqrt{(\kappa-\tilde{\kappa})^{2}+4\kappa^{+}\kappa^{-}})$
an eigenvalue of the twist matrix $K$. As a result, we can obtain constraints
on the parameters $\bar{u}$ that, up to symmetrization, provide one solution
of the Bethe equations (14). More details will be given elsewhere.

##  Appendix A Rational functions and commutations relations

Let us introduce the commutation relations of the $t_{ij}(u)$ given by

$\displaystyle[t_{ij}(u),t_{kl}(v)]=g(u,v)\big{(}t_{kj}(v)t_{il}(u)-t_{kj}(u)t_{il}(v)\big{)}.$

In particular, we will only use the following ones

$\displaystyle
t_{11}(u)t_{12}(v)=f(v,u)t_{12}(v)t_{11}(u)+g(u,v)t_{12}(u)t_{11}(v),$ (23)
$\displaystyle
t_{22}(u)t_{12}(v)=f(u,v)t_{12}(v)t_{22}(u)+g(v,u)t_{12}(u)t_{22}(v),$ (24)
$\displaystyle
t_{21}(u)t_{12}(v)=t_{12}(u)t_{21}(v)+g(u,v)\big{(}t_{11}(v)t_{22}(u)-t_{11}(u)t_{22}(v)\big{)}.$
(25)

The functions $f$ and $g$ are given by (8). They allow us to find the action
on the multiple product of $t_{12}(u)$ that forms a string of length $M$,
namely,

$\displaystyle
t_{11}(u)t_{12}(\bar{v})=f(\bar{v},u)t_{12}(\bar{v})t_{11}(u)+\sum_{i=1}^{M}g(u,v_{i})f(\bar{v}_{i},v_{i})t_{12}(u)t_{12}(\bar{v}_{i})t_{11}(v_{i}),$
$\displaystyle
t_{22}(u)t_{12}(\bar{v})=f(u,\bar{v})t_{12}(\bar{v})t_{22}(u)+\sum_{i=1}^{M}g(v_{i},u)f(v_{i},\bar{v}_{i})t_{12}(u)t_{12}(\bar{v}_{i})t_{22}(v_{i}),$
$\displaystyle t_{21}(u)t_{12}(\bar{v})=t_{12}(\bar{v})t_{21}(u)$
$\displaystyle\hphantom{t_{21}(u)t_{12}(\bar{v})=}{}+\sum_{i=1}^{M}f(u_{i},\bar{u}_{i})\big{(}g(u,u_{i})t_{22}(u_{i})t_{12}(\bar{u}_{i})t_{11}(u)+g(x_{i},y)t_{22}(u)t_{12}(\bar{u}_{i})t_{11}(u_{i})\big{)}.$

##  Appendix B Projection of the Bethe vector

The Bethe vector in terms of the operator $t_{12}(u)$ is given by

$\displaystyle\nu_{12}(\bar{u})|0\rangle=\mu^{N}\sum_{i=0}^{N}\sum_{\bar{u}\to\\{\bar{u}_{\scriptscriptstyle\rm
I},\bar{u}_{\scriptscriptstyle\rm
II}\\}}\left(\frac{\rho}{\kappa^{-}}\right)^{N-i}\,W^{N}_{i}(\bar{u}_{\scriptscriptstyle\rm
I}|\bar{u}_{\scriptscriptstyle\rm II})t_{12}(\bar{u}_{\scriptscriptstyle\rm
II})|0\rangle,$

with $\\#\bar{u}_{\scriptscriptstyle\rm II}=i$,
$\\#\bar{u}_{\scriptscriptstyle\rm I}=N-i$ a partition of the set $\bar{u}$.
The sum is over all ordered partitions, denoted
$\bar{u}\to\\{\bar{u}_{\scriptscriptstyle\rm I},\bar{u}_{\scriptscriptstyle\rm
II}\\}$. The coefficient is given by

$\displaystyle
W^{N}_{i}(u_{1},\dots,u_{N-i}|u_{N-i+1},\dots,u_{N})=\operatorname{Sym}^{N-i}_{u_{1},\dots,u_{N-i}}\left(\prod_{j=1}^{N-i}W^{N+1-j}_{N-j}(u_{j}|u_{j+1},\dots,u_{N})\right),$

where

$\displaystyle
W^{i}_{i-1}(u_{1}|u_{2},\dots,u_{i})=f(\bar{u}_{1},u_{1})\lambda_{1}(u_{1})+f(u_{1},\bar{u}_{1})\lambda_{2}(u_{1})=\Lambda^{i-1}_{d}(u_{1},\bar{u}_{1},1,1)$

and

$\displaystyle\mbox{Sym}^{M}_{\bar{u}}\big{(}F(\bar{u})\big{)}=\frac{1}{M!}\sum_{\sigma\in
S_{M}}F(\bar{u}^{\sigma}),$

with $\bar{u}^{\sigma}=\\{u_{\sigma(1)},\dots,u_{\sigma(M)}\\}$ an element of
the permutation group $S_{M}$. The dual Bethe vector in terms of the operator
$t_{21}(u)$ is given by

$\displaystyle\langle
0|\nu_{12}(\bar{u})=\mu^{N}\sum_{i=0}^{N}\sum_{\bar{u}\to\\{\bar{u}_{\scriptscriptstyle\rm
I},\bar{u}_{\scriptscriptstyle\rm
II}\\}}\left(\frac{\rho}{\kappa^{+}}\right)^{N-i}W^{N}_{i}(\bar{u}_{\scriptscriptstyle\rm
I}|\bar{u}_{\scriptscriptstyle\rm II})\langle
0|t_{21}(\bar{u}_{\scriptscriptstyle\rm II}).$

### Acknowledgements

We thank J. Avan, N. Grosjean and R Nepomechie for discussions. R.A.P. would
like to thanks the hospitality of the Laboratoire de Physique Théorique et
Modélisation at the Université de Cergy-Pontoise where a part of this work was
done. S.B. is supported by the Université de Cergy-Pontoise post doctoral
fellowship. R.A.P. is supported by Sao Paulo Research Foundation (FAPESP),
grants # 2014/00453-8 and # 2014/20364-0. We also thank the referees for their
constructive remarks.

## References

  * [1] Avan J., Belliard S., Grosjean N., Pimenta R.A., Modified algebraic Bethe ansatz for XXZ chain on the segment – III – Proof, Nuclear Phys. B 899 (2015), 229–246, arXiv:1506.0214.
  * [2] Baseilhac P., Koizumi K., Exact spectrum of the XXZ open spin chain from the $q$-Onsager algebra representation theory, J. Stat. Mech. Theory Exp. 2007 (2007), P09006, 27 pages, hep-th/0703106.
  * [3] Batchelor M.T., Baxter R.J., O’Rourke M.J., Yung C.M., Exact solution and interfacial tension of the six-vertex model with anti-periodic boundary conditions, J. Phys. A: Math. Gen. 28 (1995), 2759–2770, hep-th/9502040.
  * [4] Baxter R.J., Exactly solved models in statistical mechanics, Academic Press, Inc., London, 1982.
  * [5] Belliard S., Modified algebraic Bethe ansatz for XXZ chain on the segment – I: Triangular cases, Nuclear Phys. B 892 (2015), 1–20, arXiv:1408.4840.
  * [6] Belliard S., Crampé N., Heisenberg XXX model with general boundaries: eigenvectors from algebraic Bethe ansatz, SIGMA 9 (2013), 072, 12 pages, arXiv:1309.6165.
  * [7] Belliard S., Crampé N., Ragoucy E., Algebraic Bethe ansatz for open XXX model with triangular boundary matrices, Lett. Math. Phys. 103 (2013), 493–506, arXiv:1209.4269.
  * [8] Belliard S., Pakuliak S., Ragoucy E., Slavnov N.A., The algebraic Bethe ansatz for scalar products in ${\rm SU}(3)$-invariant integrable models, J. Stat. Mech. Theory Exp. 2012 (2012), P10017, 25 pages, arXiv:1207.0956.
  * [9] Belliard S., Pakuliak S., Ragoucy E., Slavnov N.A., Form factors in ${\rm SU}(3)$-invariant integrable models, J. Stat. Mech. Theory Exp. 2013 (2013), P04033, 16 pages, arXiv:1211.3968.
  * [10] Belliard S., Pimenta R.A., Modified algebraic Bethe ansatz for XXZ chain on the segment – II – General cases, Nuclear Phys. B 894 (2015), 527–552, arXiv:1412.7511.
  * [11] Bethe H., Zur Theorie der Metalle. I. Eigenwerte und Eigenfunktionen Atomkette, Z. Phys. 71 (1931), 205–226.
  * [12] Cao J., Yang W., Shi K., Wang Y., Off-diagonal Bethe ansatz and exact solution a topological spin rin, Phys. Rev. Lett. 111 (2013), 137201, 5 pages, arXiv:1305.7328.
  * [13] Crampé N., Algebraic Bethe ansatz for the totally asymmetric simple exclusion process with boundaries, J. Phys. A: Math. Theor. 48 (2015), 08FT01, 12 pages, arXiv:1411.7954.
  * [14] de Vega H.J., Families of commuting transfer matrices and integrable models with disorder, Nuclear Phys. B 240 (1984), 495–513.
  * [15] Derkachev S.É., The $R$-matrix factorization, $Q$-operator, and variable separation in the case of the XXX spin chain with the ${\rm SL}(2,{\mathbb{C}})$ symmetry group, Theoret. and Math. Phys. 169 (2011), 1539–1550.
  * [16] Derkachov S.É., Korchemsky G.P., Manashov A.N., Separation of variables for the quantum ${\rm SL}(2,{\mathbb{R}})$ spin chain, J. High Energy Phys. 2003 (2003), no. 7, 047, 30 pages.
  * [17] Faddeev L.D., Takhtadzhyan L.A., Spectrum and scattering of excitations in the one-dimensional isotropic Heisenberg model, J. Sov. Math. 24 (1984), 241–267.
  * [18] Galleas W., Functional relations from the Yang–Baxter algebra: eigenvalues of the XXZ model with non-diagonal twisted and open boundary conditions, Nuclear Phys. B 790 (2008), 524–542, arXiv:0708.0009.
  * [19] Gaudin M., La fonction d’onde de Bethe, Collection du Commissariat à l’Énergie Atomique: Série Scientifique, Masson, Paris, 1983.
  * [20] Gaudin M., McCoy B.M., Wu T.T., Normalization sum for the Bethe’s hypothesis wave functions of the Heisenberg–Ising chain, Phys. Rev. D 23 (1981), 417–419.
  * [21] Kitanine N., Kozlowski K.K., Maillet J.M., Niccoli G., Slavnov N.A., Terras V., Correlation functions of the open XXZ chain. I, J. Stat. Mech. Theory Exp. 2007 (2007), P10009, 37 pages, arXiv:0707.1995.
  * [22] Kitanine N., Kozlowski K.K., Maillet J.M., Niccoli G., Slavnov N.A., Terras V., Correlation functions of the open XXZ chain. II, J. Stat. Mech. Theory Exp. 2008 (2008), P07010, 33 pages, arXiv:0803.3305.
  * [23] Kitanine N., Maillet J.M., Niccoli G., Terras V., On determinant representations of scalar products and form factors in the SoV approach: the XXX case, arXiv:1506.02630.
  * [24] Kitanine N., Maillet J.M., Terras V., Form factors of the $XXZ$ Heisenberg spin-$\frac{1}{2}$ finite chain, Nuclear Phys. B 554 (1999), 647–678, math-ph/9807020.
  * [25] Korepin V.E., Calculation of norms of Bethe wave functions, Comm. Math. Phys. 86 (1982), 391–418.
  * [26] Korepin V.E., Bogoliubov N.M., Izergin A.G., Quantum inverse scattering method and correlation functions, Cambridge Monographs on Mathematical Physics, Cambridge University Press, Cambridge, 1993.
  * [27] Mukhin E., Tarasov V., Varchenko A., Bethe algebra of homogeneous XXX Heisenberg model has simple spectrum, Comm. Math. Phys. 288 (2009), 1–42, arXiv:0706.0688.
  * [28] Nepomechie R.I., Bethe ansatz solution of the open XXZ chain with nondiagonal boundary terms, J. Phys. A: Math. Gen. 37 (2004), 433–440, special issue on recent advances in the theory of quantum integrable systems, hep-th/0304092.
  * [29] Nepomechie R.I., Wang C., Algebraic Bethe ansatz for singular solutions, J. Phys. A: Math. Theor. 46 (2013), 325002, 8 pages, arXiv:1304.7978.
  * [30] Niccoli G., Terras V., Antiperiodic XXZ chains with arbitrary spins: complete eigenstate construction by functional equations in separation of variables, Lett. Math. Phys. 105 (2015), 989–1031, arXiv:1411.6488.
  * [31] Pakuliak S., Ragoucy E., Slavnov N.A., Zero modes method and form factors in quantum integrable models, Nuclear Phys. B 893 (2015), 459–481, arXiv:1412.6037.
  * [32] Ribeiro G.A.P., Martins M.J., Galleas W., Integrable ${\rm SU}(N)$ vertex models with general toroidal boundary conditions, Nuclear Phys. B 675 (2003), 567–583, nlin.SI/0308011.
  * [33] Sklyanin E.K., Quantum inverse scattering method. Selected topics, in Quantum Group and Quantum Integrable Systems, Nankai Lectures Math. Phys., World Sci. Publ., River Edge, NJ, 1992, 63–97, hep-th/9211111.
  * [34] Sklyanin E.K., Takhtadzhyan L.A., Faddeev L.D., Quantum inverse problem method. I, Theoret. and Math. Phys. 40 (1979), 688–706.
  * [35] Slavnov N.A., Calculation of scalar products of wave functions and form-factors in the framework of the algebraic Bethe ansatz, Theoret. and Math. Phys. 79 (1989), 502–508.
  * [36] Takhtadzhan L.A., Faddeev L.D., The quantum method for the inverse problem and the XYZ Heisenberg model, Russian Math. Surveys 34 (1979), no. 5, 11–68.
  * [37] Wang Y., Yang W., Cao J., Shi K., On the inhomogeneous T-Q relation for quantum integrable models, arXiv:1506.02512.


# On transaction parallelizability in Ethereum

Nadi Sarrar<EMAIL_ADDRESS>
###### Abstract

Ethereum clients execute transactions in a sequential order prescribed by the
consensus protocol. This avoids conflicts due to different transactions
reading or modifying the state of the same accounts. Sequential execution
constitutes a safe and conservative approach to blockchain transaction
processing which forgoes running transactions in parallel even when doing so
would be beneficial and safe, e.g., when there is no intersection in the sets
of accounts that the transactions touch. In particular, the possibility of
unpredictable conflicts greatly complicates the implementation of parallel
transactions. This calls for a study quantifying the occurrence of such
conflicts and analyzing the speedup potential from parallelization under such
conditions.

In this work we take a step in this direction and conduct an empirical study
based on all blocks and transactions in Ethereum to date to develop an
understanding of actual transaction execution constraints. Our results suggest
that the opportunity of safely running transactions in parallel is improving
over time. Notably, a simple scheduler already achieves average speedups of
1.77 and 1.94 in years 2017 and 2018, respectively. Using the more advanced
Heterogeneous Earliest Finish Time (HEFT) scheduler we observe even higher
gains.

We believe that there are practical ways to leverage these insights in
production, including (i) by easy parallelizability [1] and access lists [2],
and (ii) by speculatively running transactions in parallel with a mechanism to
resolve conflicts as they occur [3, 4]. The results of this work are more
closely related to (i) as we assume future knowledge of each transaction’s
access list (and processing time, in case of HEFT) and use it to determine a
schedule that never aborts.

##  1 Methodology

We use a modified version of go-ethereum
v1.8.19111https://github.com/nadisarrar/go-ethereum/ which maintains for each
transaction a record of all accounts which are relevant during its execution,
including the From account, the To account (except for create transactions),
all accounts created directly or indirectly, and all accounts that the
transaction interacts with, including value transfers, method invocations, and
any access to state such as balance and codeHash.

We feed these transaction records into three different simulators: one that is
based on transaction graph properties, which we call largest cluster, and two
transaction schedulers, simple and HEFT. We rely on gasUsed as an estimator of
transaction processing times – this is not always accurate but produces
repeatable results which are independent of specific client implementations,
hardware capabilities, network performance, etc.

###  1.1 Largest Cluster

Transactions are represented as vertices in a graph, with edges between a pair
if and only if they have at least one account in their access lists in common.
The largest cluster metric is based on the largest (in total gasUsed) disjoint
set of transactions in that graph. We calculate the total processing time as
the time it takes to sequentially process that largest disjoint subset222There
may be transactions belonging to the same transitively connected group that
can safely be processed in parallel. The largest cluster simulation therefore
leaves room for optimization.. We assume that there are infinite threads
available so that we can allocate a dedicated thread to each remaining subset
and therefore make sure that the smaller sets finish no later than the largest
one. This metric serves as a point of comparison for the other metrics that
are limited in the number of available threads.

###  1.2 Simple Scheduler

The simple scheduler processes transactions in batches of size equal to the
number of threads. Batches are allowed to contain only non-conflicting
transactions. This means that batches may contain fewer transactions than
there are threads available, and a batch only completes when its longest
running transaction completes. For these reasons some threads may be left
underutilized. Since the dataset contains access lists for each transaction,
the scheduler can determine a schedule that is guaranteed to not conflict.

###  1.3 HEFT Scheduler

HEFT [5] is a heuristics-based task scheduling algorithm that minimizes
overall completion time. The implementation of HEFT used in this study is
available at [6]. We configure HEFT to use the consensus ordering of
transactions as well as their access lists as precedence constraints, and the
gasUsed metric as an approximation of processing time. Unfortunately, the
python implementation of HEFT used in our experiments failed for some large
blocks. As a work around, we compute schedules for at most 32 transactions at
a time, which means that larger blocks require multiple runs of HEFT. We
impose the same 32 transaction limitation in the simple and largest cluster
experiments as well to ensure a fair comparison.

(a) Comparison of largest cluster and simple.

(b) Comparison with a limit of 32 transactions per schedule.

Figure 1: Parallelization gain as achieved by largest cluster, simple, and
HEFT.

##  2 Results

In Figure 1(a) we compare the performance of the largest cluster approach to
that of the simple scheduler approach. There are three main take-aways. The
direction of the LOESS regressions (solid lines) suggests that the potential
for parallelizability improves over time. However, we observe a declining
performance in 2018 which appears to coincide with a decreasing transaction
rate. We speculate that the more transactions there are, the better they can
be parallelized, due partly to there being little room for parallelization if
a block contains few transactions, and partly to an increase in the number of
popular accounts (smart contracts) and disjoint user groups generating
distinct clusters of transactions. Lastly, while the simple scheduler performs
worse than largest cluster, it still shows promising gains overall.

Figure 1(b) compares the HEFT scheduler to the other two simulations. These
experiments use a cap of 32 transactions per schedule because of a limitation
in the HEFT implementation (Sec. 1.3). This results in decreased performance
of largest cluster as well as simple. The overall trends however remain
largely the same. HEFT outperforms simple, reaching the same levels as largest
cluster.

In summary, while average thread utilization remains modest compared to the
number of available threads, all three simulations demonstrate that notable
gains are achievable with parallelization.

##  3 Future work

This study assumes that we have future knowledge of the accounts that a
transaction touches. While we may eventually be able to use access lists [2]
in an upcoming version of Ethereum for that purpose, further research on
scheduling algorithms which relax this assumption will be needed.

In case of HEFT, we further assume that gasUsed serves as an approximation of
a transaction’s processing time. This is not feasible in practice because
gasUsed only becomes available after a transaction was run, as part of the
transaction receipt. Future research may investigate the (un)predictability of
a transaction’s processing time.

This work makes no claims as to where resource bottlenecks in Ethereum
currently lie. We hope our work complements studies investigating constrained
resources such as networking and disk I/O, as well as work on alternative
concurrent state database data structures.

Finally, we hope this research direction will prove useful in the current
phases of Ethereum Serenity development, where accounts are assigned to shards
(groups), and cross-shard communication (messages between groups) incurs some
overhead. Intuitively, a greater parallelizability may suggest a lesser need
for cross-shard communication, depending on our ability to correctly identify
the clusters in real time as accounts are created and assigned to shards, and
the rate at which those assignments change.

## References

  * [1] V. Buterin, “Easy parallelizability,” Ethereum Improvement Proposal 648, June 2017.
  * [2] “Ethereum sharding specification proposal: Access lists.” https://github.com/ethereum/sharding/blob/master/docs/doc.md#access-list, accessed 01/2019.
  * [3] T. D. Dickerson, P. Gazzillo, M. Herlihy, and E. Koskinen, “Adding concurrency to smart contracts,” CoRR, vol. abs/1702.04467, 2017.
  * [4] P. S. Anjana, S. Kumari, S. Peri, S. Rathor, and A. Somani, “An efficient framework for concurrent execution of smart contracts,” CoRR, vol. abs/1809.01326, 2018.
  * [5] H. Topcuoglu, S. Hariri, and M.-Y. Wu, “Performance-effective and low-complexity task scheduling for heterogeneous computing,” IEEE Transactions on Parallel and Distributed Systems, vol. 13, March 2002.
  * [6] M. Rocklin, “HEFT implementation in Python.” https://github.com/mrocklin/heft, 2013-2017.


# Assignment of excitonic insulators in ab initio theories: the case of NiBr2

Stefano Di Sabatino<EMAIL_ADDRESS>Laboratoire de
Physique Théorique, Université de Toulouse, CNRS, UPS, F-31062 Toulouse,
France and European Theoretical Spectroscopy Facility (ETSF)  Alejandro
Molina-Sánchez<EMAIL_ADDRESS>Institute of Materials Science (ICMUV),
Universitat de València, Catedrático Beltrán 2, E-46980, Valencia, Spain  Pina
Romaniello<EMAIL_ADDRESS>Laboratoire de Physique
Théorique, Université de Toulouse, CNRS, UPS, F-31062 Toulouse, France and
European Theoretical Spectroscopy Facility (ETSF)  Davide Sangalli
<EMAIL_ADDRESS>Istituto di Struttura della Materia, Consiglio
Nazionale delle Ricerche (CNR-ISM), Division of Ultrafast Processes in
Materials (FLASHit), Via Salaria Km 29.5, CP 10, I-00016 Monterotondo
Stazione, Italy

###### Abstract

In this work we perform a detailed first–principles analysis of the electronic
and optical properties of NiBr2 within the state–of–the–art $GW$+BSE scheme to
determine whether this system displays negative excitonic energies, which
would identify it as an (half) excitonic insulator. Particular attention is
payed to the convergence of the $GW$band structure and to the consistency
between approximations employed in the ground-state calculations and
approximations employed in the linear response calculations. We show that
these two issues play a crucial role in identifying the excitonic nature of
NiBr2.

##  I Introduction

The seminal paper [1] of Jerome and coworkers predicted the realization of the
broken symmetry excitonic insulator (EI) phase in materials where the exciton
binding energy overcomes the fundamental band gap. The phase transition from a
symmetric state to the lower-symmetry EI would be purely electronically
driven. A very rich activity followed, with various research groups trying to
find signature of EI phase in different class of materials, ranging from semi-
metals to small gap semiconductors and nanostructures[2, 3, 4, 5, 6, 7, 8, 9,
10, 11, 12, 13, 14, 15].

The determination of the EI phase is mostly quantitative, rather then
qualitative, i.e. it relies on the relative value of the excitonic binding
energy versus the fundamental band gap. Therefore accurate numerical
simulations are essential. Moreover, the broken symmetry phase is often
associated to an instability of the crystal structure, and a distortion
towards a lower symmetry phase. It then becomes a chicken and egg question
whether it is the lattice distortion that drives the electronic phase
transition or vice-versa, and whether or not one should speak about excitonic
instability at all in many cases. The structural and the electronic phase
transitions cannot be disentangled experimentally, making a direct proof of
the EI existence a non trivial problem, with different experimental groups
often reaching opposite conclusions. This is why ab initio numerical
simulations are the key instrument to be used. Indeed, while indications of
the EI phase from experimental measurements have been collected, in most works
the assignment of the EI phase is based on ab initio numerical simulations,
either as support of experimental data or as purely theoretical/computational
works.

The state-of-the-art approach to model excitons in ab initio simulations is
the so-called Bethe-Salpeter equation (BSE), computed on top of an accurate
band structure including quasi-particle $GW$corrections, within the so-called
$GW$+BSE scheme. Then, numerically, a material is foreseen to be an excitonic
insulator if the BSE binding energy overcomes the $GW$band gap. This approach
has been used in several works.[7, 8, 16, 11, 10] Recently, following this
approach, a possible EI phase was identified by Jiang et al.[16] for the
magnetic layered NiBr2. NiBr2 is quite unique for two reasons. First of all,
it has a very large band gap (about 4 eV), compared to other materials
explored so far. This would suggest that the EI phase could be ruled out from
the beginning, since typical exciton binding energies range from few meV to
few hundreds meV in bulk materials. Second, NiBr2 is magnetic, thus leading to
an interplay between charge and spin degrees of freedom. As a result, the EI
phase was predicted to be realized in one spin channel only. The authors
referred to this case as an half-EI (HEI).

The $GW$+BSE scheme is a well established numerical scheme, and it works
surprisingly well in describing the optical properties of a wide range of
materials, ranging from bulk insulators and semi-conductors, to layered and 2D
materials, 1D materials, such as carbon nanotubes, up to isolated systems such
as clusters and molecules. In NiBr2, similarly to other EI candidates, the
$GW$+BSE scheme results in solutions of the excitonic hamiltonian with
slightly negative energies (few meV below zero). However the validity of the
$GW$+BSE scheme to correctly capture poles with energy close, or even below,
zero should be questioned from the beginning for various reasons. First, it is
well known that the results of the $GW$+BSE scheme depend on the starting
point, i.e. on the functional employed in the initial density functional
theory (DFT) calculation. For example, if the initial band structure is
computed starting from a DFT calculation using the PBE or the PBE0
approximation, one should refer to it as $GW$@PBE+BSE or $GW$@PBE0+BSE,
respectively. Although a change of the functional might only result in small
energy shifts, this could be enough to move an energy eigenvalue of the BSE
across zero. Second, in this protocol approximations employed in the ground-
state calculations are not consistent with the approximations used in the
linear response part. As a consequence, for example, the $GW$@DFT+BSE scheme,
if used to compute excitations in the spin flip channel, fails in reproducing
the Goldstone sum rule for magnons, i.e. it fails to give a zero energy pole
in magnetic materials [17]. This failure is independent of the DFT flavour,
the corresponding error can be large, and it clearly shows that the scheme is
not reliable for poles with energy close to zero.

Surprisingly, these issues have been discussed very little in the literature
when addressing the existence of possible EI phases. In the present manuscript
we use the quite unique case of NiBr2 to address this deficiency. We show that
NiBr2 is an extreme case where one should be very careful in proposing the
existence of an EI phase, since the state-of-the-art $GW$@DFT+BSE scheme is so
much affected by the starting point, that even qualitatively results can
change completely. These findings are not specific to the case of NiBr2 but
can be extended to the analysis of a wide range of materials.

The paper is organized as follows. In Sec. II we define in a clear way what an
EI phase (or an EI instability) is, and discuss the numerical protocols which
properly address the existence of the EI phase in general. In particular we
discuss why one has to be consistent in the approximations used to describe
the ground state and the response of the system. In Sec. III we present the
results of our numerical simulations and discuss the band structure
(subsection III.1) as well as the absorption spectrum (subsection III.2) of
ferromagnetic NiBr2 using various approximations. We finally draw our
conclusions and perspectives in Sec. IV.

##  II Theoretical protocol

The EI phase has two fundamental key properties: (i) it is a broken symmetry
phase of electronic nature, and (ii) it is described by introducing the long
range electron-hole interaction in the simulations. These are the two
essential requirements that any numerical simulation predicting an EI phase
should take into account.

Let us define the electronic hamiltonian $H$, which also includes the
potential of the atoms in a chosen crystal lattice structure, and its
eigenvalues (at fixed electron number) $E_{\lambda}$, with $E_{0}$ the ground-
state energy. We can define a corresponding point symmetry group $S$ and a
group of discrete translations $T$, such that $[H,S]=0$ and $[H,T]=0$. We call
$\Omega^{ST}$ the phase space of electronic states which respect such symmetry
groups, $\Omega$ the global phase space, and $\Omega^{S}$ ($\Omega^{T}$) the
phase space which respects the $S$ ($T$) group alone. We have
$\Omega\supseteq\Omega^{S}\supseteq\Omega^{ST}$
($\Omega\supseteq\Omega^{T}\supseteq\Omega^{ST}$). The excitonic insulator
phase arises from a spontaneous symmetry breaking, usually with respect to the
$S$ group, and will be the focus of the present manuscript. Accordingly, we
define $E^{ST}_{0}$ as the global energy minimum over $\Omega^{ST}$, and
$E^{T}_{0}$ as the global energy minimum over $\Omega^{T}$. A broken-symmetry
phase ground state exists if and only if $E^{T}_{0}<E_{0}^{ST}$.

Standard ab initio simulations compute the ground state with the symmetry
specified by the crystal lattice structure of the system, within a given
approximation for the exchange and correlation (xc) energy; we indicate the
corresponding total energy as $E_{0}^{ST,\text{xc}_{g}}$, with $\text{xc}_{g}$
referring to the xc approximation used in the ground state. There are then two
options available to check if a broken-symmetry state with energy
$E_{0}^{T,\text{xc}_{g}}<E_{0}^{ST,\text{xc}_{g}}$ exists: (i.1) to perform a
second ground-state simulation without symmetry or (i.2) to compute the
spectrum of the neutral excitations of the material,
$\omega_{\lambda}=(E_{\lambda}-E_{0})$, at zero transferred momentum.111The
finite momentum excitations would describe the breaking of the transnational
symmetries as well. We do not consider this case here. In the latter case
there is no need to impose symmetry breaking. Indeed, neutral excitations are
computed from the solution of the linear response function $\chi(\omega)$. The
excited-state energies, $E_{\lambda}$, and the corresponding wavefunctions,
$\Psi_{\lambda}$, which appear in the Lehmann representation of
$\chi(\omega)$, do not belong to $\Omega_{ST}$ in general. More precisely,
given the space symmetry group $S$, the states $\Psi_{\lambda}$ can be
classified as degenerate multiplets of the irreducible representations
$R_{S}^{\lambda}$ of the space group. Symmetry operations of $\Omega_{ST}$ can
in general send a state of a given multiplet into any other state of the same
multiplet. Because of that, the states of a monodimensional representation
$R_{S}^{\lambda}$ are singlets, i.e. nondegenerate, and belong to
$\Omega_{ST}$. Instead, states of a multidimensional representations
$R_{S}^{I}$ do not belong to $\Omega_{ST}$, and they must be considered
together to obtain a closed sub-space with respect to the symmetry operation
of the group $S$. If a multiplet of states has energy lower than the ground
state, then a spontaneous symmetry breaking would be favored, with an emerging
ground state composed by a random linear combination of the states in the
multiplet.

If option (i.1) is chosen, the same level of approximation must be used in the
ground-state calculation with and without symmetry. Comparing
$E_{0}^{T,\text{xc}_{g1}}$ and $E_{0}^{ST,\text{xc}_{g2}}$, with
$\text{xc}_{g1}$ and $\text{xc}_{g2}$ two different xc approximations used in
the ground state, does not have any validity. Moreover, the long range
electron-hole interaction needs to be considered. In some works the EI phase
has been predicted on the basis of calculations using the local density
approximation (LDA) or the generalized gradient approximation (GGA), i.e.
$\text{xc}_{g}=\text{LDA/GGA}$. While also within these approximations a
symmetry-broken phase could be obtained, such phases should not be labelled as
EI, since LDA/GGA completely miss the long-range behaviour of the
electron–hole interaction. Indeed, it is well known that time dependent DFT
(TDDFT), within the LDA or GGA approximations (TDLDA or TDGGA), cannot
describe excitons in the excited states, and thus these approximations cannot
capture the EI phase in the ground state either.

In case option (i.2) is chosen, the total energy of possible broken-symmetry
phases can be obtained by defining
$E_{\lambda}^{T,\text{xc}}=E_{0}^{ST,\text{xc}_{g}}+\omega^{ST,\text{xc}_{e}}_{\lambda}$,
with
$\omega^{ST,\text{xc}}_{\lambda}=E_{\lambda}^{T,\text{xc}}-E_{0}^{ST,\text{xc}}$
a zero-momentum neutral excitation of the system, using consistent
approximations between ground state ($\text{xc}_{g}$) and excited state
($\text{xc}_{e}$) calculations. Then the criteria to establish the existence,
for some $\lambda_{0}$, of a broken-symmetry phase becomes
$\omega^{ST,\text{xc}_{e}}_{\lambda_{0}}<0$. In such situation we can identify
$E_{\lambda_{0}}^{T,\text{xc}}=E_{0}^{T,\text{xc}}$ as the global ground state
energy. The state-of-the-art scheme to describe excitons in the ab initio
framework is BSE. $\omega^{\text{xc}_{e}}_{\lambda}<0$, with ${\text{xc}}_{e}$
referring to the xc approximation used for the BSE kernel, is then often used
as the golden standard to predict the existence of the excitonic insulator
phase. Within this approach the essential consistency condition at the core of
option (i.2) would rely on the use of a one-body Green’s function which is a
self-consistent solution of the Dyson equation for a given self-energy
approximation and on the use of the same level of approximation for the self-
energy and the BSE kernel [17].

However, in practical calculations, this consistency condition is violated
since the state-of-the-art is the $GW$@DFT+BSE scheme. Indeed, the BSE kernel
is obtained from the Hartree plus screened-exchange self-energy (HSEX)
approximation, i.e. $\text{xc}_{e}=\text{HSEX}$, while
$\text{xc}_{g}=GW\text{@DFT}$. The self-consistent DFT ground state with $GW$
corrections should be replaced by a self-consistent simulation with
$\text{xc}_{g}=\text{HSEX}$. Even ignoring this point and assuming that $GW$
is a good starting point for a BSE calculation, at least an eigenvalue self-
consistent $GW$ (ev$GW$) procedure should be carried out. However both SC-HSEX
and ev$GW$ can be very demanding. An alternative path is to use DFT with
hybrid functionals, which are computationally affordable. Hybrid functional
simulations include (a fraction of) the long–range electron–hole interaction
and can be used to capture EI phases. In this case
$\text{xc}_{g}=\text{xc}_{e}=\text{Hybrid}$ can be used, by employing TDDFT
for the excited states. To underline this consistency, we will use here the
notation PBE0+TDPBE0 for example for the PBE0 case. Indeed, it has been shown
that hybrid functionals can capture excitonic excited states [18] and EI
phases [9]. In general the exchange fraction entering hybrids is not properly
screened, and this could be particularly problematic in 2D materials, since
the 2D nature of the screening is completely missing. Nevertheless, as we
shall see, it remains a valid approach to demonstrate the concept of
consistency we discussed above.

In this work we will adopt the protocol described in option (i.2) to identify
possible EI phases in NiBr2 and we will calculate the optical spectrum within
$GW$@PBE+BSE, $GW$@PBE0+BSE and PBE0+TDPBE0. This will allow us to address
both the starting-point dependence of the $GW$@DFT+BSE scheme and to consider
a proper scheme with $\text{xc}_{g}=\text{xc}_{e}$.

##  III Results

We studied NiBr2 in its ferromagnetic ground state. Density functional theory
calculations were performed using the Quantum Espresso code [19] within both
the Perdew-Burke-Ernzerhof (PBE) and the PBE0 exchange correlation functional
[20]. Optimized norm-conserving Vanderbilt pseudopotentials have been
generated using ONCVPSP [21]. The self-consistent (scf) simulations converged
with an energy cutoff of 80 Ry and a $\mathbf{k}$-point grid of $9\times
9\times 1$ (the direct band gap at $\Gamma$ has been converged to 0.01 eV).

We then used the DFT wave-functions $\psi^{KS}_{n\mathbf{k}}(\mathbf{r})$ and
energies $\epsilon^{KS}_{n\mathbf{k}}$ to perform $GW$ and BSE simulations.
For the PBE case we re-computed $\psi^{KS}_{n\mathbf{k}}(\mathbf{r})$ and
$\epsilon^{KS}_{n\mathbf{k}}$ with a non self-consistent (nscf) calculation on
a $19\times 19\times 1$ $\mathbf{k}$-point grid in the Brillouin zone. For the
PBE0 case instead we used a $9\times 9\times 1$ $\mathbf{k}$-point grid. This
is due both to the higher computational cost of PBE0, and to the fact that,
with hybrid functionals, it is not possible, in the nscf calculations, to use
a $\mathbf{k}$-point grid different from the one used for the scf case. As we
shall see, this restriction does not affect our conclusions. We computed the
one-shot $GW$quasiparticle (QP) corrections on top of DFT using the Yambo code
[22, 23]. Dynamical screening effects have been taken into account using the
Godby-Needs plasmon-pole approximation (PPA) [24]. Both in DFT and,
successively, in $GW$ and BSE, we used a slab model with a 40 a.u. vacuum
thickness. Moreover, the Coulomb cut-off technique [25] has been used to
ensure no interactions between periodic images. In the $GW$step, we used a
cutoff of 40 Ry for the wave-functions; an energy cut-off four times bigger,
i.e. of 320 Ry, for $V_{\text{xc}}$ and for the exchange part of the Self-
Energy $\Sigma_{x}$; finally an energy cutoff of 10 Ry for the screening and
the correlation parts of the self-energy $\Sigma_{c}$ (see the Supplemental
Material more for details [26]). We included a total of 200 bands in the
screening and 500 in the evaluation of the self-energy. In the optical
response we took into account excitonic effects via the Bethe-Salpeter
Equation [27, 28]. The BSE absorption spectrum has been calculated using a
static screening converged using the same parameters of the $GW$calculation.
We used a cutoff of 4 Ry for the electron-hole exchange and 10 Ry for the
screened interaction. We included bands from 7 to 23 to obtain a converged
spectrum up to 10 eV.

In the following we demonstrate that the occurrence of negative energy
solutions of the BSE strongly depends on the choice of the approximations used
in the ground state and for the BSE kernel. In particular the lack of self-
consistency in the $GW$ band structure and of consistency between ground-state
approximations and BSE kernel approximations can lead to a questionable EI
phase.

Figure 1: Spin-resolved band structure (majority- and minority-spin channels
are reported in blue and red, respectively) for $GW$@PBE (left) and $GW$@PBE0
(right) are reported with a continuous line. The DFT band structures used as
starting point for the $GW$correction are reported with a dashed line. The VBM
is set to zero.

###  III.1 Band structure

The starting point for a BSE calculation is a single-particle band structure
obtained correcting the DFT energies. The approximation chosen for the DFT
determines the quality of the screened interaction $W$, which is computed
within the RPA approximation using Kohn-Sham energies. This results in the so-
called $W_{0}$ screened interaction. $W_{0}$ enters the one-shot $GW$or
$G_{0}W_{0}$ approach, which is usually employed for the band structure. The
one-shot $GW$is also employed in the present work. Even without invoking self-
consistency, $GW$calculations can be problematic in 2D materials. We found,
indeed, that, for the system under study, particular attention has to be paid
when converging the $GW$ quasiparticle corrections. See also the discussion in
the Supplemental Material [26] (see, also, references [29, 30] therein).

| PBE | PBE0 | $GW$
---|---|---|---
|  | [30] |  | @PBE | @PBE [30] | @PBE0
ind. gap M-m | 0.83 | 0.88 | 3.79 | 2.82 | 2.62 | 5.75
dir. gap m-m | 1.33 | 1.34 | 4.45 | 2.91 | 2.62 | 6.56
dir. gap M-M | 4.09 | 3.95 | 6.52 | 5.01 | 5.37 | 7.41
dir. gap M-m | 0.99 | 1.04 | 3.86 | 3.05 | 3.06 | 5.91
Table 1: Indirect and direct band gaps (in eV).

In Fig. 1 we show PBE and PBE0 spin-resolved band structures, separated in
majority (M) and minority (m) spin channels. The corresponding $GW$@PBE and
$GW$@PBE0 band structures are also shown. In all cases the band gap is
indirect, with 0.83 eV at the PBE level compared to 3.79 eV at the PBE0 level,
and it is due to transitions from the majority valence band maximum (VBM) at
$\Gamma$ towards the minority conduction band minimum (CBM), which is at
roughly half way along the $\Gamma$M high symmetry line. Also the direct band
gaps within PBE0 (4.45 eV and 6.52 eV for m-m and M-M spin channels,
respectively) are significantly larger than the ones obtained using PBE (1.33
eV and 4.09 eV for m-m and M-M spin channels, respectively). This is expected
due to the fraction of HF exchange taken into account in the PBE0 functional.

Moving to $GW$, one would naively expect that, despite the starting-point
dependence, the $GW$@PBE and $GW$@PBE0 gaps should be close to each other,
with $GW$@PBE approaching the expected self-consistent $GW$ result from below,
i.e. with the $GW$ corrections opening the PBE gap, and with $GW$@PBE0 from
above, i.e. with the $GW$ corrections closing the PBE0 gap. However this is
clearly not the case as shown in Fig. 1. $GW$@PBE opens the PBE gap from 0.83
to 2.62 eV. Direct band gaps are increased to 2.91 eV and 5.01 eV for the
minority and majority spin channels, respectively. However $GW$@PBE0 also
opens the PBE0 gap from 3.79 to 5.75 eV. Direct band gaps are 6.56 eV and 7.41
eV for the minority and majority spin channels respectively. See also Tab. 1
for a comparison of all values. This clearly points to the fact that
$W_{0}^{\text{PBE}}$ and $W_{0}^{\text{PBE0}}$ are different and that,
accordingly, the two $G_{0}W_{0}$ self-energies are different. One can then
wonder which result is the most reliable.

On one hand the $GW$@PBE is the state-of-the-art scheme to compute quasi-
particle band structure. Indeed, the PBE screening entering
$W_{0}^{\text{PBE}}$ is usually a reasonable approximation to the physical
screening due to error cancellations between (i) the underestimated
fundamental gap, and (ii) the neglected excitonic effects. On the other hand
hybrid functionals like PBE0 are expected to better approximate quasi-particle
wave-functions and energies. Also simulations based on the hybrid HSE0
functional [31, 16] give a band structure rather close to the PBE0 band-
structure. However $GW$@hybrids is a scheme which has not been much tested
yet, since the compatibility of $GW$ codes with hybrid functionals is a quite
recent feature, and, despite being very promising [32, 33, 34], they are not
very well tested. By correcting the electronic gap, hybrid functionals tend to
underestimate the electronic screening which may result in a overestimation of
the $GW$@hybrids gap [32]. Indeed, in NiBr2 the underestimated screening leads
to a very strong exchange term, even bigger than the exchange fraction
entering in PBE0. As a consequence, it is more reasonable to directly compare
$GW$@PBE with hybrids like PBE0. The estimated experimental band gap, probed
using STS on NiBr2 thin films, down to the monolayer, ranges between 3.4 and
4.5 eV, in good agreement with the PBE0 value [35].

Beside this evident, but yet quantitative difference, there is a more
important qualitative difference between PBE and PBE0 band structures, which
is also inherited by the corresponding $GW$ band structure. Within PBE there
exist three almost flat bands in the minority channel which are very close in
energy to the top of the valence band in the majority spin channel. These
bands are Nickel $3d$ states [16], in particular with $t_{2g}$ character,
which are likely not well described within PBE. Indeed, they are not present
within our PBE0 simulations, nor in the HSE0 and LDA+$U$ results in the
literature [31, 16]. Moving to $GW$@PBE these same states remain there,
pointing to the fact that the PBE wave-functions are not well described, and
that PBE is not a good starting point. These states also show a quite peculiar
behaviour in the $GW$ simulations, and, as we shall see, they are responsible
for the very sharp difference in the spectrum of $GW$@PBE+BSE compared to
PBE0+TDPBE0 or $GW$@PBE0+BSE. In the Supplemental Material [26] we discuss the
convergence of the $GW$@PBE quasi-particle corrections. We highlight the
“unusual” behaviour of these flat Ni $t_{2g\downarrow}$ bands, and we discuss
in detail the very slow convergence compared to the energy cut-off used in the
description of the screened interaction. The discussion agrees with previous
numerical studies [36], although this fact has been somehow overlooked in past
simulations on 2D materials [16]. Indeed $GW$@PBE corrections on these
occupied Ni $t_{2g\downarrow}$ states lead to a situation in which the
majority and the minority VBM are very close in energy (during the
$GW$convergence their relative position is even inverted). The direct gap
switches character, from M-m to m-m, and the system almost becomes a direct
gap material with an energy difference between $GW$@PBE direct and indirect
gap of 9 meV (compared to 500 meV in PBE, 660 meV in PBE0 and 810 meV in
$GW$@PBE0). The small m-m gap is also part of the reason why negative poles
will appear in the response function calculated at the level of $GW$@PBE+BSE,
as we will discuss in Sec. III.2.

The conclusion of this analysis is that NiBr2 is an indirect band gap
semiconductor. The indirect gap is majority-minority, as expected in standard
magnetic materials. Looking at the direct gap, this is also majority-minority,
the minority-minority gap is slightly larger (less than 1 eV), while the
majority-majority gap is significantly larger. This behaviour is due to the
empty spin-minority Ni $3d$ orbitals, with $e_{g}$ character, which lay in
between the more predominantly $p$-like valence and conduction bands.

###  III.2 Optical properties

Figure 2: Imaginary part of the macroscopic polarizability calculated within
$GW$@PBE+BSE (upper panel), $GW$@PBE0+BSE (middle panel), and PBE0+TDPBE0
(lower panel). The results obtained within the corresponding independent
particle approximation, i.e. ignoring electron-hole interaction, are also
reported (dashed lines). Curves are vertically shifted for clarity. Bright and
dark excitons energies are reported as vertical thicks.

We now turn our attention to the optical properties of NiBr2. In Fig. 2 the
$GW$@PBE+BSE absorption spectrum is reported together with the $GW$@PBE0+BSE
and the PBE0+TDPBE0 spectra. All spectra are dominated by an intense excitonic
peak, which we label $E_{0}$. All methods describe $E_{0}$ as an exciton
mostly located nearby $\Gamma$ in the BZ (see Fig. 3) and corresponding to a
transition from $p$ bands $\epsilon_{p\mathbf{k}\downarrow}$ to $d(e_{g})$
bands $\epsilon_{d\mathbf{k}\downarrow}$ in the spin minority channel. At
energy larger then $E_{0}$ the spectral intensity is strongly reduced and
eventually it goes up again at even larger energy. The position of the $E_{0}$
peak shifts from 3.3 eV ($GW$@PBE+BSE), to 3.6 eV (PBE0+TDPBE0), and to 4.1 eV
($GW$@PBE0+BSE). The $E_{0}$ energy is always lower compared to the
$(\epsilon_{d\mathbf{k}\downarrow}-\epsilon_{p\mathbf{k}\downarrow})$ energy
difference. We define its binding energy
$b_{E_{0}}=\min_{\mathbf{k}}(\epsilon_{d\mathbf{k}\downarrow}-\epsilon_{p\mathbf{k}\downarrow})-E_{0}$.
$b_{E_{0}}\approx 1$ eV in both $GW$@PBE+BSE (1.1 eV) and PBE0+TDPBE0 (0.9
eV), while $b_{E_{0}}=2.4$ eV for $GW$@PBE0+BSE. This discrepancy is again due
to the underestimated RPA screening obtained using the PBE0 starting point.
This leads to an overestimation of the band-gap which is here compensated by
an overestimation of the binding energy and results in a peak position not
very different from $GW$@PBE+BSE and PBE0+TDPBE0. As for the band structure
case, the conclusion here is that it is better to directly compare
$GW$@PBE+BSE to PBE0+TDPBE0.

The qualitative behaviour of the spectrum is instead significantly different
at energy below $E_{0}$. All methods show poles at energy well below $E_{0}$,
due to transitions between occupied and flat bands in the minority spin
channel, which we identify as $d(t_{2g})$ states (see Supplemental Material
[26]), and the empty spin minority $d(e_{g})$-states. However, within
$GW$@PBE+BSE, this leads to a series of three bright excitons to appear, which
we label $D_{i}$ in Fig. 2 following the choice done in Ref. [16]. These peaks
are located at 1.67 eV ($D_{1}$), 1.98 eV ($D_{2}$), 2.21 eV ($D_{3}$). As
already pointed out, the description of the occupied $t_{2g}$ states is
questionable in PBE, and so are the absorption features related to them.
Indeed the series of the $D$ peaks is missing within $GW$@PBE0+BSE and
PBE0+TD-BPE0, where all $d(t_{2g})\rightarrow d(e_{g})$ transitions are dark.
This is probably due to the fact that the PBE starting point mixes $d$ and $p$
bands character, while PBE0 does not.

Figure 3: Band contribution to the main excitonic peak $E_{0}$: (left panel)
$GW$@PBE+BSE, (middle panel) $GW$@PBE0+BSE, (right panel) PBE0+TDPBE0. Note
that the plotted band structures are the DFT ones. We also report the
excitonic wavefunction in reciprocal space.

Beside these aspects, we are here interested in discussing whether or not
NiBr2 is an HEI. We thus need to turn our attention to the lowest energy
solutions among the $d(t_{2g})\rightarrow d(e_{g})$ transitions. As discussed,
the occupied $d(t_{2g})$ states are very close to the valence band maximum
within $GW$@PBE, while they are deeper in energy within PBE0 and $GW$@PBE0. As
a consequence the $GW$@PBE+BSE approach produces negative eigenvalues. There
is a doubly degenerate dark exciton at -0.17 eV ($X_{1,1}$) and a non-
degenerate exciton at -0.16 eV ($X_{1,2}$). These two excitons share a similar
physics and are labeled as $X_{1}$ in Fig. 2. Other two dark excitons located
at 0.59 ($X_{2,1}$) and 0.68 eV ($X_{2,2}$) are indicated with $X_{2}$. Both
$X_{1}$ and $X_{2}$ are associated with transitions between $t_{2g}$ and
$e_{g}$ bands in the minority channel (see Supplemental Material [26]). The
$X_{1}$ poles in particular are the ones which were used to claim the possible
existence of an HEI ground state in Ref. [16]. Instead, within $GW$@PBE0+BSE
and PBE0+TD-BPE0 no negative eigenvalues appear and therefore there is no EI
phase. Since the qualitative behaviour of these three approaches is the same,
and, as discussed in the introduction, the existence of a possible EI phase is
mostly quantitative, we have now to address the question: which result is the
most reliable? Both $GW$@DFT+BSE schemes are based on a one-shot
$GW$calculation and plagued by inconsistencies between the approximation used
for the ground state and the approximation used for the excited states. This
points to the fact that, as far as this quantitative question is concerned,
PBE0+TDPBE0 is the scheme to rely on. For NiBr2 it should be noticed that the
$GW$@DFT+BSE inconsistency is not so drastic. Indeed the physics of both the
ground state and the excited states is determined by the same screened
interaction $W$. However, the state-of-the-art $GW$@PBE+BSE is affected by a
bad description of the occupied $d$ bands in the current study, which is here
crucial to address this quantitative answer. $GW$@PBE0+BSE suffers from a PBE0
screening which is likely underestimated, although errors induced in the band
structure are cancelled by the overestimation of the binding energy, thus
affecting in a less severe way the final position of the BSE peaks. We also
notice that in Ref. [16], the HEI phase is also supported by HSE0+BSE
simulations where, however, the inconsistency between ground state and excited
states approximation is more problematic.

In conclusion, our analysis of the band structure and the of the optical
spectral features indicates that PBE0+TDPBE0 gives a good overall description
of NiBr2. Since this scheme does not produce any negative-energy pole in the
optical absorption, we conclude that NiBr2 is not an HEI. This is the most
important result of this work.

##  IV Conclusions and perspectives

In this work we addressed the study of excitonic insulator phases within
common ab initio theoretical approaches such as time-dependent density
functional theory and the Bethe-Salpeter equation of Many-Body Perturbation
Theory. We discussed the fact that in order to correctly capture this phase a
theoretical approach should be able to describe two important features, i.e.
an electronic-driven broken symmetry and the long range electron-hole
interaction. That is why TDDFT within the common LDA/GGA approximations (which
miss the long range electron-hole interaction) cannot correctly describe this
phase. The BSE is more appropriate, but, in order to correctly describe an
electronic-driven broken symmetry, it should be used within consistent
approximations. We discussed that this corresponds to use the same level of
approximation for the self-energy used in a self-consistent ground-state
calculation and for the kernel of the BSE.

For the specific case of ferromagnetic NiBr2, by comparing the optical
spectrum and the energy position of the lowest energy dark excitons obtained
with different schemes, we arrived at the conclusion that the PBE0+TDPBE0
scheme is the most reliable, also because it uses the same level of
approximation for the ground state and for the excited states. Within
PBE0+TDPBE0 all dark poles have a positive energy, which indicates that NiBr2
is not an HEI, contrary to previous assignments. Moreover, we find that the
occurrence of negative energy peaks in the $GW$@PBE+BSE is due to the starting
point dependence of the $GW$ scheme, rather than to the inconsistency between
ground-state and excited state approximations. Indeed, the PBE starting point
gives a questionable description of the $t_{2g}$ states, which leads, within
$GW$@PBE, to $t_{2g}$ bands located too high with respect to the occupied $p$
states. This does not happen within PBE0 and $GW$@PBE0 (and also with respect
to other approximations studied in the literature such as HSE0 and PBE+$U$).

More in general, our study points to the fact that using consistent
approximations in ab initio theories is of crucial importance in order to have
reliable results.

###### Acknowledgements.

S. Di Sabatino and P. Romaniello would like to acknowledge financial support
by the EUR grant NanoX ANR-17-EURE-0009 in the framework of the “Programme des
Investissements d’Avenir” and by the ANR (project ANR-18-CE30-0025 and
ANR-19-CE30-0011). A. M.-S. acknowledge financial support by Ramón y Cajal
programme (grant RYC2018-024024-I; MINECO, Spain), Agencia Estatal de
Investigación (AEI), through the project PID2020-112507GB-I00 (Novel quantum
states in heterostructures of 2D materials), and Generalitat Valenciana,
program SEJIGENT (reference 2021/034), project Magnons in magnetic 2D
materials for a novel electronics (2D MAGNONICS), and Planes complementarios
de I+D+I en materiales avanzados, project SPINO2D, reference MFA/2022/009.

## References

  * Jérome  _et al._ [1967] D. Jérome, T. M. Rice, and W. Kohn, Phys. Rev. 158, 462 (1967).
  * Bucher  _et al._ [1991] B. Bucher, P. Steiner, and P. Wachter, Phys. Rev. Lett. 67, 2717 (1991).
  * Bronold and Fehske [2006] F. X. Bronold and H. Fehske, Phys. Rev. B 74, 165107 (2006).
  * Cercellier  _et al._ [2007] H. Cercellier, C. Monney, F. Clerc, C. Battaglia, L. Despont, M. G. Garnier, H. Beck, P. Aebi, L. Patthey, H. Berger, and L. Forró, Phys. Rev. Lett. 99, 146403 (2007).
  * Su and MacDonald [2008] J.-J. Su and A. H. MacDonald, Nature Physics 4, 799 (2008).
  * Wakisaka  _et al._ [2009] Y. Wakisaka, T. Sudayama, K. Takubo, T. Mizokawa, M. Arita, H. Namatame, M. Taniguchi, N. Katayama, M. Nohara, and H. Takagi, Phys. Rev. Lett. 103, 026402 (2009).
  * Ataei  _et al._ [2021] S. S. Ataei, D. Varsano, E. Molinari, and M. Rontani, Proceedings of the National Academy of Sciences 118, 10.1073/pnas.2010110118 (2021), https://www.pnas.org/content/118/13/e2010110118.full.pdf .
  * Varsano  _et al._ [2017] D. Varsano, S. Sorella, D. Sangalli, M. Barborini, S. Corni, E. Molinari, and M. Rontani, Nature Communications 8, 1461 (2017).
  * Hellgren  _et al._ [2018] M. Hellgren, J. Baima, and A. Acheche, Phys. Rev. B 98, 201103 (2018).
  * Dong and Li [2020] S. Dong and Y. Li, Phys. Rev. B 102, 155119 (2020).
  * Jiang  _et al._ [2020] Z. Jiang, W. Lou, Y. Liu, Y. Li, H. Song, K. Chang, W. Duan, and S. Zhang, Phys. Rev. Lett. 124, 166401 (2020).
  * Mazza  _et al._ [2020] G. Mazza, M. Rösner, L. Windgätter, S. Latini, H. Hübener, A. J. Millis, A. Rubio, and A. Georges, Phys. Rev. Lett. 124, 197601 (2020).
  * Baldini  _et al._ [2020] E. Baldini, A. Zong, D. Choi, C. Lee, M. H. Michael, L. Windgaetter, I. I. Mazin, S. Latini, D. Azoury, B. Lv, A. Kogar, Y. Wang, Y. Lu, T. Takayama, H. Takagi, A. J. Millis, A. Rubio, E. Demler, and N. Gedik, The spontaneous symmetry breaking in ta2nise5 is structural in nature (2020).
  * Sun  _et al._ [2022] B. Sun, W. Zhao, T. Palomaki, Z. Fei, E. Runburg, P. Malinowski, X. Huang, J. Cenker, Y.-T. Cui, J.-H. Chu, X. Xu, S. S. Ataei, D. Varsano, M. Palummo, E. Molinari, M. Rontani, and D. H. Cobden, Nature Physics 18, 94 (2022).
  * Chen  _et al._ [2022] D. Chen, Z. Lian, X. Huang, Y. Su, M. Rashetnia, L. Ma, L. Yan, M. Blei, L. Xiang, T. Taniguchi, K. Watanabe, S. Tongay, D. Smirnov, Z. Wang, C. Zhang, Y.-T. Cui, and S.-F. Shi, Nature Physics 18, 1171 (2022).
  * Jiang  _et al._ [2019] Z. Jiang, Y. Li, W. Duan, and S. Zhang, Phys. Rev. Lett. 122, 236402 (2019).
  * Müller [2017] M. C. T. D. Müller,  _Spin-wave excitations and electron-magnon scattering in elementary ferromagnets from ab initio many-body perturbation theory_ , Ph.D. thesis (2017).
  * Paier  _et al._ [2008] J. Paier, M. Marsman, and G. Kresse, Phys. Rev. B 78, 121201 (2008).
  * Giannozzi  _et al._ [2009] P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, G. L. Chiarotti, M. Cococcioni, I. Dabo, A. D. Corso, S. de Gironcoli, S. Fabris, G. Fratesi, R. Gebauer, U. Gerstmann, C. Gougoussis, A. Kokalj, M. Lazzeri, L. Martin-Samos, N. Marzari, F. Mauri, R. Mazzarello, S. Paolini, A. Pasquarello, L. Paulatto, C. Sbraccia, S. Scandolo, G. Sclauzero, A. P. Seitsonen, A. Smogunov, P. Umari, and R. M. Wentzcovitch, Journal of Physics: Condensed Matter 21, 395502 (2009).
  * Perdew  _et al._ [1996] J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996).
  * Hamann [2013] D. R. Hamann, Phys. Rev. B 88, 085117 (2013).
  * Marini  _et al._ [2009] A. Marini, C. Hogan, M. Grüning, and D. Varsano, Computer Physics Communications 180, 1392 (2009).
  * Sangalli  _et al._ [2019] D. Sangalli, A. Ferretti, H. Miranda, C. Attaccalite, I. Marri, E. Cannuccia, P. Melo, M. Marsili, F. Paleari, A. Marrazzo, G. Prandini, P. Bonfà, M. O. Atambo, F. Affinito, M. Palummo, A. Molina-Sánchez, C. Hogan, M. Grüning, D. Varsano, and A. Marini, Journal of Physics: Condensed Matter 31, 325902 (2019).
  * Godby and Needs [1989] R. W. Godby and R. J. Needs, Phys. Rev. Lett. 62, 1169 (1989).
  * Rozzi  _et al._ [2006] C. A. Rozzi, D. Varsano, A. Marini, E. K. U. Gross, and A. Rubio, Phys. Rev. B 73, 205119 (2006).
  * [26] See Supplemental Material at [URL will be inserted by publisher].
  * Martin  _et al._ [2016] R. M. Martin, L. Reining, and D. M. Ceperley, Interacting electrons: Theory and computational approaches (Cambridge University Press, 2016).
  * Rohlfing and Louie [2000] M. Rohlfing and S. G. Louie, Phys. Rev. B 62, 4927 (2000).
  * Guandalini  _et al._ [2022] A. Guandalini, P. D’Amico, A. Ferretti, and D. Varsano, arXiv preprint arXiv:2205.11946 (2022).
  * [30] Computational 2D materials database (C2DB).
  * Mushtaq  _et al._ [2017] M. Mushtaq, Y. Zhou, and X. Xiang, RSC Adv. 7, 22541 (2017).
  * Fuchs  _et al._ [2007] F. Fuchs, J. Furthmüller, F. Bechstedt, M. Shishkin, and G. Kresse, Phys. Rev. B 76, 115109 (2007).
  * Leppert  _et al._ [2019] L. Leppert, T. Rangel, and J. B. Neaton, Phys. Rev. Mater. 3, 103803 (2019).
  * Gant  _et al._ [2022] S. E. Gant, J. B. Haber, M. R. Filip, F. Sagredo, D. Wing, G. Ohad, L. Kronik, and J. B. Neaton, Phys. Rev. Mater. 6, 053802 (2022).
  * Bikaljević  _et al._ [2021] D. Bikaljević, C. González-Orellana, M. Peña-Diaz, D. Steiner, J. Dreiser, P. Gargiani, M. Foerster, M. A. Niño, L. Aballe, S. Ruiz-Gomez, N. Friedrich, J. Hieulle, L. Jingcheng, M. Ilyn, C. Rogero, and J. I. Pascual, ACS Nano 15, 14985 (2021), pMID: 34491033, https://doi.org/10.1021/acsnano.1c05221 .
  * Hüser  _et al._ [2013] F. Hüser, T. Olsen, and K. S. Thygesen, Phys. Rev. B 88, 245309 (2013).