# HubbardNet: Efficient Predictions of the Bose-Hubbard Model Spectrum with
Deep Neural Networks

Ziyan Zhu  Present address: Stanford Institute for Materials and Energy
Sciences, SLAC National Accelerator Laboratory, Menlo Park, CA 94025, USA
<EMAIL_ADDRESS>Department of Physics, Harvard University, Cambridge,
Massachusetts 02138, USA  Marios Mattheakis  John A. Paulson School of
Engineering and Applied Sciences, Harvard University, Cambridge, Massachusetts
02138, USA  Weiwei Pan  John A. Paulson School of Engineering and Applied
Sciences, Harvard University, Cambridge, Massachusetts 02138, USA  Efthimios
Kaxiras  Department of Physics, Harvard University, Cambridge, Massachusetts
02138, USA  John A. Paulson School of Engineering and Applied Sciences,
Harvard University, Cambridge, Massachusetts 02138, USA

###### Abstract

We present a deep neural network (DNN)-based model (HubbardNet) to
variationally find the ground state and excited state wavefunctions of the
one-dimensional and two-dimensional Bose-Hubbard model. Using this model for a
square lattice with $M$ sites, we obtain the energy spectrum as an analytical
function of the on-site Coulomb repulsion, $U$, and the total number of
particles, $N$, from a single training. This approach bypasses the need to
solve a new hamiltonian for each different set of values $(U,N)$. Using
HubbardNet, we identify the two ground state phases of the Bose-Hubbard model
(Mott insulator and superfluid). We show that the DNN-parametrized solutions
are in excellent agreement with results from the exact diagonalization of the
hamiltonian, and it outperforms exact diagonalization in terms of
computational scaling. These advantages suggest that our model is promising
for efficient and accurate computation of exact phase diagrams of many-body
lattice hamiltonians.

##  I Introduction

In condensed matter physics, simple many-body lattice models have been
employed to capture the essential physics of intriguing phenomena like phase
transitions, strongly correlated states, and exotic phases. A particularly
popular and powerful simple model is the Hubbard model, which has proven
useful in elucidating, for instance, the phase diagram and the spectra of
high-temperature superconductors [1, 2, 3], of unconventional
superconductivity in twisted multilayered van der Waals heterostructures [4,
5, 6, 7, 8, 9, 10, 11, 12], and the superfluid-to-insulator transition in
bosonic gases [13, 14, 15, 16]. However, the solution of the deceptively
simple Hubbard model remains a challenge due to the fact that, as a function
of the system size and the number of particles in the system, an exponentially
large space is required to fully describe the many-body wavefunction.

In this work, we adopt a variational approach by considering a class of
wavefunctions parameterized by a deep neural network (DNN) to reduce the
parameter space to search for the optimal solution. Similar machine learning
(ML)-based variational approaches have proven useful in providing insights
into the ground state of interacting spin system [17, 18, 19], bosonic systems
[20, 21, 22], fermionic systems [23, 24, 25, 26, 24, 25], dynamics in many-
body systems [27], and open quantum systems [28, 29]. While the amount of
information encoded in these ML-based models in principle scales linearly as a
function of the number of particles and sites, the actual training of a DNN is
known to be computationally costly and in many cases more costly than the
exact diagonalization (ED) of the hamiltonian, because it requires a new
training for every different value of the parameters on which the hamiltonian
depends.

The very large size of the space needed to represent the wavefunction has so
far restricted the search for solutions, typically, only to the ground state
of the many-body system. Our proposed model (in the following referred to as
the HubbardNet) employs a variational approach to solve for the ground state
as well as for the excited states of the one-dimensional (1D) and two-
dimensional (2D) Bose-Hubbard model. Our DNN-parametrized solutions are an
analytic function of the on-site Coulomb repulsion between electrons,
described by the parameter $U$, and the total number of particles $N$. This
bypasses the need to perform an optimization for every value of the $(U,N)$
parameters of the hamiltonian. Our results provide the proof of principle that
variational approach employing DNN’s can be used to efficiently map out the
phase boundaries in $(U,N)$ space of low-dimensional lattice models.

The paper is organized as follows. In Section II we introduce the Bose-Hubbard
model and its construction. In Section III we present the details of
HubbardNet and use it to solve for the ground state and excited states of the
1D and 2D Bose-Hubbard model on a square lattice with $M$ sites, using both
periodic and open boundary conditions. We evaluate the performance of our
method and compare to ED in Section IV. Finally, we summarize our findings in
Section V and discuss the possible extensions to fermionic systems.

##  II The Bose-Hubbard Model

The Bose-Hubbard Hamiltonian is given by

$\hat{\mathcal{H}}=-t\sum_{\langle
ij\rangle}\hat{a}_{i}\hat{a}_{j}^{\dagger}+\frac{U}{2}\sum_{i=1}^{M}\hat{n}_{i}(\hat{n}_{i}-1)+\sum_{i}^{M}(V_{i}-\mu)n_{i},$
(1)

where the sum on $i$ runs over the $M$ lattice sites and $\langle ij\rangle$
denotes the nearest neighbor pairs in the lattice. The first term in the
hamiltonian represents the kinetic energy of the quantum particles, with
$\hat{a}_{i}^{\dagger},\hat{a}_{i}$ being the creation and annihilation
operators for a particle at the lattice site $i$, and $-t$ being the amount of
energy gained by a particle during a “hop” between two neighboring sites
$\langle ij\rangle$. The next term in the hamiltonian represents the
interaction between particles, assumed to be charged (electrons or ions). We
assume the system consists of only one kind of charged particles, with $U$
describing the Coulomb repulsion between each pair;
$\hat{n}=\hat{a}_{i}^{\dagger}\hat{a}_{i}$ is the number operator, whose
expectation value gives the occupation of the lattice site $i$. The last term
includes a contribution from an external potential $V_{i}$ at site $i$, and
the chemical potential energy of the system, $\mu$, which we will assume to be
zero in this work. In the following we take $t=1$ and express all values of
energy in units of $t$. This simple model can capture the superfluid phase for
$U/t\leq 1$ and the Mott insulator phase for $U/t\gg 1$ [11].

We will denote a state, referred to as a Fock state, of the system as
$|\bm{n}_{j}\rangle$, which is defined by the occupation of each site
$n_{i}^{(j)}$, with the index $i$ running over all $M$ sites, that is,

$|\bm{n}_{j}\rangle=|n_{1}^{(j)},n_{2}^{(j)},\dots,n_{M}^{(j)}\rangle,$ (2)

with $n_{i}^{(j)}\in[0,N]$ for $i=1,\dots,M$, which satisfies

$\sum_{i=1}^{M}n_{i}^{(j)}=N.$ (3)

The number $\mathcal{N}_{B}$ of these Fock states for the Bose-Hubbard model
with $N$ particles in $M$ sites is given by

$\mathcal{N}_{B}={M+N-1\choose N}.$ (4)

Using this basis of Fock states, the many-body wavefunction $|\Psi_{k}\rangle$
that describes the entire system is written as:

$|\Psi_{k}\rangle=\sum_{j=1}^{\mathcal{N}_{B}}\psi_{j}^{(k)}|\bm{n}_{j}\rangle,$
(5)

and the solution involves finding all the many-body wavefunctions that satisfy
the eigenvalue equation

$\hat{\mathcal{H}}|\Psi_{k}\rangle=E_{k}|\Psi_{k}\rangle,$ (6)

with $E_{k}$ being the $k$-th energy eigenvalues. The states
$|\Psi_{k}\rangle$ form an orthonormal set,
$\langle\Psi_{i}|\Psi_{j}\rangle=\delta_{ij}$, which implies

$\sum_{j=1}^{\mathcal{N}_{B}}|\psi_{j}^{(k)}|^{2}=1$ (7)

for each value of $k$. Thus, the solution involves finding the complex numbers
$\psi_{j}^{(k)},j=1,\dots,\mathcal{N}_{B}$, that describe each of the many-
body states.

The solution is found by constructing the hamiltonian matrix

$\mathcal{H}_{ij}=\langle\bm{n}_{i}|\hat{\mathcal{H}}|\bm{n}_{j}\rangle,\;\;i,j=1,\dots,\mathcal{N}_{B},$
(8)

and finding all its eigenvalues and eigenvectors. To this end, every Fock
state needs to be associated with an integer label, in a unique and known
ordering. This can be done iteratively by Ponomarev ordering [30, 31], a
procedure described in detail by Reventós et al. [32]. We note that the
kinetic-energy term produces contributions only to the off-diagonal
hamiltonian matrix elements, because it involves hops of particles between
nearest-neighbor sites, which changes the occupations of these sites. For
example, for a model with $M=5,N=3$, supposing the $i^{\mathrm{th}}$ state is
$\bm{n}_{i}=[2,1,0,0,0]$ and the $j^{\mathrm{th}}$ state is
$\bm{n}_{j}=[1,2,0,0,0]$ by Ponomarev ordering, $\mathcal{H}_{ij}=-t$. This is
because a particle at site 1 in state $\bm{n}_{i}$ can hop to site 2, which
results in state $\bm{n}_{j}$. This nearest neighbor hop contributes a term of
$-t$ to the energy. By similar arguments, the Coulomb repulsion term produces
contributions to the diagonal matrix elements only, because it does not
involve any changes in the site occupations.

Figure 1:  Schematic diagram of the data-free neural network to solve for the
Bose-Hubbard model. The network takes $(M+2)$ inputs, which are the number of
particles $n_{i}$ at all $i=1,\dots,M$ sites as well as $U$, the Coulomb
repulsion parameter, and $N$, the total number of particles in the system (in
the example shown we use $M=4$ for illustration purposes). The output layer is
two-dimensional, with the outputs $u_{1}$ and $u_{2}$ corresponding to the
real and imaginary part of the wavefunction component corresponding to this
particular configuration of particles.

##  III Neural Network-Parametrized Solution

We parametrize the many-body wavefunction $|\Psi_{k}\rangle$ by a fully
connected DNN

$|\Psi_{k}\rangle=\sum_{j=1}^{\mathcal{N}_{B}}g(\bm{n}_{j},U,N;\bm{W}^{(k)})|\bm{n}_{j}\rangle$
(9)

where $\bm{W}{(k)}$ denotes the network parameters for this state; a schematic
representation of the network is shown in Fig. 1. Thus, the network takes as
input the pair of values $(U,N)$ and the set of occupation numbers
$n_{i}^{(j)},i=1,\dots,M$, that correspond to a specific Fock state
$|\bm{n}_{j}\rangle$, and gives as output the coefficient $\psi_{j}^{(k)}$
that corresponds to this Fock state, which is a component of the many-body
wavefunction $|\Psi_{k}\rangle$ with energy $E_{k}$. There is apparently a
redundancy in these definitions, because of the constraint given in Eq. (3),
so not all of the input values to the DNN are independent. However, a better
description is that for a given value of $N$, only the Fock states within the
space defined by Eqs. (2) and (3) are allowed as inputs to the network. In
this sense, $N$, which is an important variable in the problem, can be
considered as an independent variable. To obtain the full wavefunction
$|\Psi_{k}\rangle$, the trained network needs to be applied to each of the
Fock states once.

The size of the hamiltonian matrix, $\mathcal{N}_{B}\times\mathcal{N}_{B}$,
increases as the factorials of the number of particles $N$ and sites $M$. The
exact diagonalization of the hamiltonian matrix quickly becomes intractable as
the system size increases. In contrast, the number of input parameters scales
linearly with $M$ in the DNN approach [17]. Another important advantage of the
DNN-parameterized solution is that with $U$ and $N$ being the input parameters
to the network, the network becomes an analytic function of $U$ and $N$ and we
essentially solve for a class of hamiltonians during a single training. The
network can then be used to make predictions at different values of $U$ that
are not in the training set. In contrast, other techniques including ED and
determinant quantum Monte Carlo [33] require a new calculation for every
different set of values $(U,N)$.

In this work, we use 4 hidden layers and use a hyperbolic tangent, $\tanh()$,
activation function for all hidden layers. In the examples that we considered
in this work, the wavefunctions are real. However, the wavefunctions can in
general be complex, and thus we need two real-number outputs for the real and
imaginary parts of the coefficients

$\psi_{j}^{(k)}=u_{1,j}^{(k)}+iu_{2,j}^{(k)}.$ (10)

For the ground state wavefunction, $|\Psi_{0}\rangle$ (denoted by the index
$0$), we choose the activation function of the output layer to be an
exponential, namely

$\psi_{j}^{(0)}=\exp\left(u_{1,j}^{(0)}+iu_{2,j}^{(0)}\right).$ (11)

This choice is made for better convergence with uniformly sampled initial
weights. Note that with the exponential function, the wavefunction
$\psi_{j}(\bm{n})$ is positive. The coefficients of the Fock states have the
same sign for the ground state, which is guaranteed by the structure of the
hamiltonian for $-t<0$ on all the off-diagonal matrix elements. For excited
states, we take a linear activation function for the output layer, given in
Eq. (10), because the exponential activation function would no longer work, as
the wavefunction needs to be orthogonal to the ground state wavefunction, and
the coefficients $\psi_{j}^{(k)}$ must have mixed signs.

The average of any operator $\hat{\mathcal{A}}$ that applies to Fock states
$|\bm{n}_{j}\rangle$ can be obtained through the expression

$\displaystyle\langle\hat{\mathcal{A}}\rangle_{k}\equiv\langle\Psi_{k}|\hat{\mathcal{A}}|\Psi_{k}\rangle=\sum_{i,j=1}^{\mathcal{N}_{B}}\psi_{i}^{(k)*}\langle\bm{n}_{i}|\hat{\mathcal{A}}|\bm{n}_{j}\rangle\psi_{j}^{(k)}.$
(12)

Using $\hat{\mathcal{H}}$ as the operator in Eq. (12) we can obtain the energy
$E_{k}$ of the eigenstate $|\Psi_{k}\rangle$ of the hamiltonian. In the
network representation of the many-body wavefunction, Eq. (9), the energy
takes the form:

$\displaystyle E_{k}(U,N)=$
$\displaystyle\sum_{i,j=1}^{\mathcal{N}_{B}}g^{*}(\bm{n}_{i},U,N;\bm{W}^{(k)})\mathcal{H}_{ij}g(\bm{n}_{j},U,N;\bm{W}^{(k)})$
$\displaystyle\times\left[\sum_{i=1}^{\mathcal{N}_{B}}|g(\bm{n}_{i},U,N;\bm{W}^{(k)})|^{2}\right]^{-1},$
(13)

where $\mathcal{H}_{ij}$ are the hamiltonian matrix elements defined in Eq.
(8). In the above equation we have indicated explicitly the dependence of the
energy eigenvalues on $U$ and $N$. Accordingly, we choose the loss function as
follows:

$\mathcal{L}_{k}=\frac{1}{\mathcal{N}_{U}\mathcal{N}_{N}}\sum_{l=1}^{\mathcal{N}_{U}}\sum_{m=1}^{\mathcal{N}_{N}}E_{k}(U_{l},N_{m}),$
(14)

where $E_{k}(U_{l},N_{m})$ is the expression defined in Eq. (13), and
$\mathcal{N}_{U}$, $\mathcal{N}_{N}$ are the numbers of values of $U$ and $N$
included in the training, respectively. The loss functions for the ground-
state energy $E_{0}$ and for excited-state energies $E_{k}(k>0)$, need to be
treated differently, as described in detail next.

Figure 2: One-dimensional Bose-Hubbard model solution with a periodic boundary
condition and $M=6,N=5$ obtained using the HubbardNet. (a) Trace of the loss
function offset by the negative of the minimum total energy. (b) Comparison
between the ground truth ground state energies obtained from ED (black solid
line) and the neural network results for training set (red scattered crosses)
and prediction (blue scattered points). (c)-(d) Absolute value of the
wavefunction $|\psi_{0}(\bm{n})|$ from ED (black solid lines) and neural
network (red scattered crosses) for (c) $U=5.0$ (d) $U=7.5$. (e)-(f) The
occupation number $\langle n_{i}\rangle_{0}$ at different sites that
correspond to (c)-(d). Red solid lines represent a constant occupation number
$\langle n_{i}\rangle=N/M=0.833$.

###  III.1 Ground state

To obtain the ground state energy and wavefunction, we minimize the loss as
defined in Eq. (14) directly. As an example, we first discuss the solution for
a 1D system with periodic boundary conditions. In this case, every site has 2
nearest neighbors. We train the network with three different $U$ values,
namely $U=2.0,5.0,8.0$, for a system with $M=6,N=5$, which yields
$\mathcal{N}_{B}=252$ Fock states. Here and in the rest of the paper, unless
otherwise specified, we choose the network width to be $D=400$. We use
stochastic gradient descent (SGD) for the optimization of the network weights
$\bm{W}^{(0)}$ that minimize the loss function. We use momentum for learning
$=0.9$. We adopt a cosine annealing scheme for the learning rate [34], with
the maximum learning rate chosen to be $0.01$ and decaying to zero. We reset
the learning rate every 1000 iterations. We train the network until the
variance of the loss function during the last 200 steps is less than $1\times
10^{-7}$ (see Fig. 2(a)).

In Fig. 2(b), we compare the ground-state energy values as a function of $U$,
as obtained by HubbardNet with the values obtained by diagonalizing of the
hamiltonian matrix (namely ED). The comparison includes the values of $U$ that
are part of the training set (red scattered points), as well as values beyond
the training set (blue scattered crosses). We note that even though the energy
is a smooth function of $U$, the output of the network includes all the
wavefunction components $\psi_{j}^{(0)},j=1,\dots,\mathcal{N}_{B}$, which are
a more complex function of the parameter $U$. For example, in Fig. 2(c)-(e) we
compare the wavefunction components (labeled with an arbitrary ordering to
reflect the symmetry of the basis) from the neural network output, shown as a
scatter plot, and the corresponding values obtained from ED, for three
different values of $U$, only one of which was included in the training set,
namely $U=5.0$ (see Fig. 2(d)). In all three cases, there is an excellent
agreement between the two results, indicating that with as little as
$\mathcal{N}_{U}=3$ points in the training sets, HubbardNet is capable of
predicting the wavefunctions and energies of the 1D Bose-Hubbard model
accurately for a wide range of the values of $U$ In Fig. 2(f), we show the
average occupation number $\langle\hat{n}_{i}\rangle_{0}$, obtained from Eq.
(12), using $\hat{n}_{i}$ as the operator, with the trained weights, for three
different values of $U$, both within and beyond the training set. As expected
from the periodic boundary conditions, and the large value of $U$, the average
occupation number is uniform for all values of $U$ in this example. In other
words, the large repulsion between the charged particles uniformly distributes
them in the lattice.

Figure 3: Two-dimensional Bose-Hubbard model solution with an open boundary
condition and $M=16,N=3$ with ($M_{x}=4,M_{y}=4$) obtained using the
HubbardNet. (a) Trace of the loss function offset by the negative of the
minimum total energy. (b) Comparison between the ground truth ground state
energies obtained from exact diagonalization (black solid line) and the neural
network results for the training set (red scattered points) and prediction
(blue scattered points). (c) Absolute value of the wavefunction
$|\psi_{0}(\bm{n})|$ from the exact diagonalization (solid black lines) and
neural network (red scattered crosses) for $U=2.5$. (d) The average particle
occupation number $\langle n_{i}\rangle_{0}$ corresponding to (c) with
$U=2.5$.

We next proceed to examine the capability of HubbardNet to handle more
complicated cases, namely lattices with a larger number of sites and
particles. We consider a 2D Hubbard model with open boundary conditions (a
cluster of lattice sites), with $M=16$ (a $4\times 4$ square lattice) and
$N=3$. These choices lead to a large number of Fock states,
$\mathcal{N}_{B}=816$. The hamiltonian construction is similar to the 1D case,
but now every site has 4 nearest neighbors instead of 2, except for the sites
on the boundary which have only 2 or 3 neighbors, from the open boundary
condition. In this way, the hamiltonian becomes denser compared to the 1D
case.

In Fig. 3 we show the solution, using the same values of $U$ as the training
set, namely $U=2.0,3.5,5.0$. Despite the increased size of the problem,
HubbardNet performs extremely well: both the energies, shown in Fig. 3(b), and
the wavefunction components, shown in Fig. 3(c), match the results from the ED
solution for $U=2.5$, despite the fact that this value was not part of the
training set. The occupation number distribution (Fig. 3(d)) is no longer
uniform but is instead concentrated in the center of the cluster due to the
open boundary condition.

Figure 4: One-dimensional Bose-Hubbard model ground state solution with a
periodic boundary condition for $M=5$, and the training set including both the
$N=5$ and $N=4$ states. (a) Trace of the loss function offset by the negative
of the minimum total energy. (b) Energy vs. $U$ from ED (solid black lines)
and HubbardNet (scattered points) with $N=4$ (red) and $N=5$ (blue)
respectively. Black circles are from the training set. (c)-(d) Comparison
between wavefunction components $|\psi_{0}(\bm{n})|$ obtained from ED (blue
solid lines) and HubbardNet (black solid lines) for (c) $N=4$ and (d) $N=5$
with $U=2.5$.

In addition to including several values of $U$ as part of the network input,
the number of particles in the system, $N$, can also be part of the network
input. Although the space of Fock states has different size for different
values of $N$, HubbardNet can still be used as long as $M$ is fixed, since the
energy is an average over all particle configurations in the Fock basis. In
Fig. 4, we show the results from using HubbardNet for a system with multiple
values of both $N$ and $U$ as inputs, so that the output is a function of both
variables. Figure 4(b) shows that the two curves describing the energy as a
function of $U$ for the two different values of $N$ cross each other, making
the the optimization much more difficult. Remarkably, energies obtained from
HubbardNet both at the values of $U$ included in the training set and outside
of the training set agree with the ED results (scattered points in Fig. 4),
despite the crossover between the $N=4$ and $N=5$ energy vs. $U$ curve. The
HubbardNet wavefunction components, while retaining the approximate shape of
the ED result, are slightly different from the latter (Fig. 4(c)-(d)). It is
more difficult for the network to train with multiple values of $N$ in the
input because each $N$ gives rise to a different Hilbert space size, and thus
the corresponding wavefunctions are expected to be more different than in the
case with only multiple values of $U$. Another difficulty is that $N$ needs to
be an integer, but we are not constraining the network to predict integer
values of $N$. Therefore, the parameter space for the minimization is
unnecessarily large, which makes it more difficult for the network to
converge.

###  III.2 Identifying the ground state phases

The Bose-Hubbard model with on-site Coulomb and nearest-neighbor coupling
interactions has two phases: the superfluid phase and the Mott insulator
phase. The phase transition is exhibited by a commensurate number of particles
$N/M\in\mathbb{N}$. For the 1D model, the phase transition is of the
Berezinskii-Kosterlitz-Thouless (BKT) type [35], which is very sensitive to
the finite system size, making it difficult to pinpoint the phase boundary in
a finite-size system. Since HubbardNet enables the efficient prediction of the
wavefunctions and energy values for different $U/t$ and $N$, we can use it to
identify the existence of the two phases and approximate the phase boundary.
In this section, we first apply a site-dependent external potential to mimic
an optical lattice and study the distribution of average particle number
$\langle n_{i}\rangle_{0}$, and we then examine the finite-size scaling of the
insulating gap for the 1D Bose-Hubbard model.

Figure 5: Energy and average particle number obtained with HubbardNet for
$M=N=7$ with a periodic boundary condition and a site-dependent potential with
the strength $V=1$, trained with $U=1.0,4.0,7.0,10.0,13.0$. (a) Comparison
between the ground state energy HubbardNet (scattered points) and the ED
(solid line). Red scattered points are in the training set and blue scattered
points are predictions. (b) The average particle number per site $\langle
n_{i}\rangle_{0}$ for different $t/U$ values. Black dashed lines are from the
training set.

We choose a site-dependent external potential of the form

$V_{j}=V(j-(M-1)/2)^{2}\;\;forj=0,\dots,M-1$ (15)

with $V=1$, which mimics the weak harmonic potential in ultracold atom
experiments imposed due to the profile of the laser beam [36, 20]. Fig. 5(a)
shows the results of training with 5 different values of $U$ across a wide
range; the energies from the training set and the testing set agree well with
ED. Fig. 5(b) shows the particle redistribution for different values of the
$t/U$ ratio. At a small $U$, the particles redistribute according to the
applied potential with a peak in the center of the system. At a large $U$,
they distribute uniformly as the strong Coulomb repulsion prohibits them from
freely hopping, which is the signature of a Mott insulator. We show in Fig.
5(b) that with as little as five points in the training set, we are able to
accurately predict this qualitative trend for a wide range of $U$ values.

The upper and lower boundaries of the Mott-insulator lobe are defined as the
excitation energy of one particle and one hole, $\mu^{\pm}$ respectively. For
filling factor $N/M=1$ and a fixed system size $M$, $\mu^{\pm}$ are given as,

$\displaystyle\mu^{+}(t/U)=E_{0}(t/U,M+1)-E_{0}(t/U,M),$
$\displaystyle\mu^{-}(t/U)=E_{0}(t/U,M)-E_{0}(t/U,M-1).$ (16)

The critical point of the Mott-superfluid transition is the value of $t/U$ at
which the upper and lower boundaries of the Mott lobe coincide with each other
in the thermodynamic limit [35].

In Fig. 6(a), we plot $\mu^{\pm}$ as a function of $t/U$ for $M=4,\dots,7$,
which shows the Mott lobe and is similar to the ED result obtained in Raventós
_et al._ [32]. For each $M$, we perform a single training with $N=M-1,M,M+1$
and $U=4.0,7.0,10.0$. For a finite system size, the values of $\mu^{\pm}$ are
not expected to coincide [32]. Therefore, we calculate $\mu^{\pm}$ as a
function of $M$ and note that $\mu^{\pm}$ approach each other as $M$ gets
larger. At a sufficiently large $M$, we expect $\mu^{\pm}$ to asymptotically
approach the black line in Fig. 6, which is obtained with dynamical density
matrix renormalization group theory with $M=128$ [37].

The single-particle gap of the Mott insulator phase is defined as the
difference between $\mu^{\pm}$:

$\displaystyle\Delta_{M}(t/U)$ $\displaystyle=\mu^{+}(t/U)-\mu^{-}(t/U)$ (17)

Fig. 6(b) shows $\Delta_{M}/U$ for the same values of $M$ as in 6(a); for
increasing $t/U$, the gap size decreases as expected from the Mott to
superfluid transition. In an infinite system, we expect the gap to close for
the superfluid state, with the $t/U$ value at which the gap closes being the
critical point. We can estimate the critical point value $(t/U)_{c}$ following
a similar fitting procedure as in Kashurnikov  _et al._ [38], namely, by
performing a polynomial fit to $\Delta_{M}/U$ as a function of $1/M$ for each
value of $t/U$ and finding the $t/U$ value at which $\Delta_{M}/U$ stops
varying. Kashurnikov  _et al._ [38] and [32] fit $\Delta_{M}/U$ with a degree
5 polynomial. Here, as a crude estimate, we perform a linear fit to
$\Delta_{M}/U$ as a function of $1/M$ at each $t/U$ due to the limited number
of values of $M$ at our disposal. In this way, $\Delta_{\infty}/U$ is the
intercept of the fit. From this procedure, we estimate the critical value to
be $(t/U)_{c}\approx 0.320\pm 0.013$. In comparison, the value obtained by the
determinant quantum Monte Carlo (DQMC) method is $(t/U)_{c}=0.275\pm 0.005$
[38] and the value obtained with the DMRG method is $0.305\pm 0.001$ [39].
While we our value does not agree with these other two values within the error
bars, our result is only meant to provide a crude estimate with a rather
limited number of data points.

Figure 6: Mott lobe and particle-hole excitation gap obtained with HubbardNet.
(a) The upper ($\mu^{+}$) and lower ($\mu^{-}$) boundaries of the Mott lobe
for $M=N=4,\dots,8$ obtained with HubbardNet. For each $M$, we train with
three different values of $N$ with $N=M-1,M,M+1$ and three different values of
$U$ with $U=4.0,7.0,10.0$. Black solid lines are obtained using dynamical
density-matrix renormalization group theory (DDMRG) with $M=128$, adapted from
Ejima  _et al._ [37]. (b) Charge gap in the unit of $U$ obtained with
HubbardNet for the same values of $M$ as (a).

###  III.3 Excited states

In order to obtain excited states with the approach described above, we employ
the fact that the eigenstates of the hamiltonian are orthogonal to each other,
that is, $\langle\Psi_{i}|\Psi_{j}\rangle=\delta_{ij}$. An excited state is a
local minimum of the energy. Therefore, we can obtain excited states by
minimizing the expectation value in Eq. (13) for state $|\Psi_{k}\rangle$ with
the constraint that it is orthogonal to
$|\Psi_{0}\rangle,...,|\Psi_{k-1}\rangle$. The constraint is achieved
iteratively through the Gram-Schmidt orthogonalization process. To facilitate
the calculations, we express the wavefunctions $|\Psi_{k}\rangle$ as vectors
$\bm{v}_{k}$ of dimension $\mathcal{N}_{B}$, with entries $\psi_{j}^{(k)}$. In
this notation, assume that $\bm{v}_{0},\bm{v}_{1},...,\bm{v}_{k-1}$ are the
vectors that represent the output directly from the neural network, which are
linearly independent, with $|\Psi_{0}\rangle=\bm{v}_{0}$ being the ground
state wavefunction, and the rest being the wavefunctions of excited states.
The excited state wavefunctions are obtained iteratively as follows:

$|\Psi_{k}\rangle=\bm{v}_{k}-\sum_{j=0}^{k-1}\frac{\bm{v}_{j}\cdot\bm{v}_{k}}{\bm{v}_{j}\cdot\bm{v}_{j}}\bm{v}_{j}.$
(18)

This expression for state $\bm{v}_{k}$ incorporates the projection out of this
state of components that correspond to all previous states, $j=0,\dots,k-1$.
The projection is performed at every iteration and the set of $\bm{v}_{j}$
vectors is used to compute the loss function. In other words, for an excited
state $k$, we train a new network with the outputs being constrained to be
orthogonal to the space of all states with lower energy, from the ground state
to the $(k-1)$ excited state. With the excited state optimization for multiple
values of $U$ as the input, we add the following penalty term to the loss
function:

$\displaystyle\mathcal{L}_{P}=-\exp\left(-|\bar{E}_{k}-\bar{E}_{0}|\right),$
(19)

where the index $k$ labels the excited state, $\bar{E}_{k}$ is the average
energy over the values of $U$ as obtained at each epoch and $\bar{E}_{0}$ is
the averaged ground state energy. This choice of the loss function encourages
the excited-state energy $E_{k}$ to be as close to the ground state energy as
possible. Similarly to the ground state, we can train the network as a
function of $U$ for a given excited state $k$. In Fig. 7(a) we show the energy
spectra for 4 different values of $U$, namely $U=2.0,3.0,4.0,5.0$, with one
training for a given $k$, for $M=5,N=3$ with an open boundary condition.

Figure 7: HubbardNet energy spectrum trained using $U=2.0,3.0,4.0,5.0$ as the
input, with L2-regularization with the coefficient $=1\times 10^{-4}$ and
perturbing the $U$ by an amplitude $0.01$. (a) Energy from the HubbardNet
(scattered points) vs. ED (dashed lines) for the training set (b) Energy
spectrum for different values of $U$ in the testing set (scatter points).
Eigenvalues from ED are in black solid lines.  Figure 8: Comparison between
excited state wavefunctions from ED and HubbardNet for the $4^{\mathrm{th}}$
(top row), the $15^{\mathrm{th}}$ (middle row), and the $30^{\mathrm{th}}$
(bottom row) excited states for $U=1.5$ (left) and $U=5.0$ (right). In all
subpanels, black dashed lines are from ED and red scattered points are from
HubbardNet.

We apply an L2 regularization by adding the squared sum of the network weights
to the loss function with a coefficient $10^{-4}$ in order to discourage
energies from being drastically different for similar values of $U$, an issue
that we observed in the absence of such regularization. In addition, since our
training set only contains a very small number of values of $U$ but we would
like to predict energies and wavefunctions as a smooth function of all values
of $U$, we randomly perturb the value of the input $U$ during every epoch,
with a perturbation amplitude of 0.01. In Fig. 7(a) we show that for the $U$
values in the training set, there is an excellent agreement between the
results obtained with HubbardNet and those obtained with ED. Using the
network, we can plot the full spectrum of the Bose-Hubbard model for the given
$M$ and $N$ and an arbitrary value of $U$. In Fig. 7(b) we show the spectrum
produced using the trained network at inference. The spectrum qualitatively
agrees with the ED results, except for $U\lesssim 2$ and at high energies.
This is expected because $U<2$ corresponds to values beyond the range of
training-set values, and higher-energy states are more difficult to obtain
from the network due to error accumulation in the iterative orthogonalization
scheme, Eq. (18).

We note that there is a gap in the energy spectrum at approximately $U=2.0$,
separating a few of the high-energy excited states from the rest. It is also
possible that the network has more difficulty in making accurate predictions
near the phase boundary, which by definition means that a small perturbation
in the value of $U$ leads to a large change in the qualitative behavior, and
thus a denser sampling is more desirable near such values of $U$.

In Fig. 8, we show the corresponding wavefunctions for three different excited
states. For the results presented in this figure, only $U=5.0$ is included in
the training set. For excited states with relatively low energy, for example,
the 4th excited state, the predicted wavefunctions agree with the ground truth
even for values of $U$ not in the test set. However, for excited states with
higher energy, such as the 15th or the 30th state, there is a significant
difference between the HubbardNet-predicted wavefunctions and the ground
truth, especially for $U\lesssim 2$. The discrepancy in the wavefunction is
also reflected in the corresponding energy spectrum in Fig. 7.

As an alternative approach to the state-by-state iterative optimization, one
could treat the state $k$ as an input to the network and minimize the total
energy of multiple states simultaneously. With this approach, we first project
the neural network onto the nullspace of the ground state and minimize the
average energy over $k$ states in the nullspace. In this case, the loss
function landscape becomes very complex and has many degenerate local minima.
While this approach avoids the error accumulation of the iterative
optimization approach, we found that the network produces less accurate
predictions and the results depend more sensitively on the training
parameters.

##  IV Performance Evaluation

In all the simulations presented in Sec. III, we used a single Nvidia K80 GPU
and most of the calculations can be carried out on CPUs, which suggests that
our approach is relatively computationally inexpensive. In Fig. 9, we compare
the performance of obtaining the ground state energy for a single value of $U$
using ED and the HubbardNet method, performed on a 2.7-GHz Quad-Core Intel
Core i7. Here, we take $M=N$ and a network width $D=200$. We adopt an
optimizer constant learning rate $=0.01$ without L2 regularization, and our
convergence criterion is the variance of the loss smaller than $1\times
10^{-6}$ for all cases. In Fig. 9(a), we show that the ground state energies
obtained from ED and HubbardNet are in excellent agreement for all cases. In
Fig. 9(b) we show the computational time as a function of the system size:
HubbardNet has a better scaling compared to the ED method, even though the
latter method is faster for smaller systems, with crossover being in the range
$\mathcal{N}_{B}\sim 5000$. Since the goal of numerical studies of the Hubbard
model is to reach as large systems as possible, the advantage of the
HubbardNet method is obvious.

An alternative approach to further improve the scaling of the DNN-based
approach is that instead of calculating the sum over every pair of
$\bm{n}_{i},\bm{n}_{j}$ vectors in Eq. (12), one could evaluate the
expectation value using a Monte Carlo sampler, which has been implemented in
the package netket [40]. However, due to the stochastic nature of the Monte
Carlo sample, we observed that the energy to deviate from the ground truth.
The accuracy of predictions is expected to deteriorate even more for excited
states of higher energy due to the error accumulation in the iterative scheme.

Figure 9: Performance comparison between ED and HubbardNet using a 1D chain
(for $M=N=3,\dots,8$) and periodic boundary conditions, as a function of the
system basis-size $\mathcal{N}_{B}$, see Eq. (4). (a) Percentage error of the
ground state energies obtained with HubbardNet versus ED (ground truth). (b)
Computation time for HubbardNet (black) and ED (red).

##  V Discussion and Outlook

In addition to the scaling advantage discussed in the previous section, our
approach provides the analytical solution for a family of hamiltonians with
different $U$ and $N$ values without the need to diagonalize a new hamiltonian
with every change of each parameter. This suggests that after one training
over a set of values of $U$ and $N$, the network could in principle predict
the phase diagram efficiently. We used HubbardNet to show that it is capable
of identifying the ground state phases and provides the correct finite-size
scaling for the Mott insulating gap (Sec. III.2 ). The ability to accurately
locate the phase transition point is still limited by the training time, which
can be improved by implementing more efficient architectures such as a
convolutional neural networks [21] and by using a Monte Carlo sampler [17, 20,
40]. We provide a proof-of-principle demonstration in Fig. 7(b) that after
training on a very small training set, we were able to obtain the energy
spectrum of the system as a function of $U$. In particular, the network is
able to predict the gap opening near $U=2.0$. As we point out in Sec. III.3,
energies and wavefunctions from HubbardNet deviate at the gap opening point,
and we suggest that this problem can be solved by increasing the training set
size. The choice of training set size is a tradeoff between accuracy and
training time: the larger the training set, the more accurate the predictions
should be, but the harder it is to converge and the more effort is needed to
tune the values of the hyperparameters. We also note that while the scaling of
HubbardNet is promising, the performance evaluation we showed is only for a
single state. The task of obtaining all excited states is significantly more
computationally costly due to the iterative nature of our approach.
Nevertheless, for most practical purposes, it is sufficient to compute only
the lowest few energy states, in which case HubbardNet would be an appropriate
choice. Obtaining the first few excited states allows us to compute the finite
temperature properties and identify the critical temperature. For example, for
a given temperature $T$, we can compute an observable by the Boltzmann
distribution:

$\hat{\mathcal{O}}(\beta)=\frac{1}{\mathcal{Z}}\sum_{k}\langle\Psi_{k}|\hat{\mathcal{O}}|\Psi_{k}\rangle
e^{-\beta(E_{k}-E_{0})},$ (20)

where $\beta=1/k_{B}T$, $k_{B}$ is the Boltzmann constant, $E_{k}$ is the
energy of the $k$-th excite state, and
$\mathcal{Z}=\sum_{k}\exp[-\beta(E_{k}-E_{0})]$ is the partition function.
Computing the temperature dependence of observables and extracting the
critical temperatures from HubbardNet merits future work.

While we only provided here a proof-of-concept demonstration of HubbardNet,
this approach could be easily adapted to solve more realistic problems
including longer-ranged hopping parametrized according to the specific
materials’ properties as well as systems in the presence of an external
potential. In addition, in this work, we solve the Bose-Hubbard model in real
space and obtained the wavefunction components in the Fock basis. We could
perform a Bloch expansion of Eq. (1) in momentum space and variationally
minimize the expectation of the hamiltonian, with the momentum $\bm{k}$ taken
as an input to the network, which enables very fast prediction of the band
structure.

An interesting extension of the HubbardNet approach would be to use it for the
study of fermionic systems and for obtaining the phase diagrams of low-
dimensional quantum systems including twisted bilayer van der Waals
heterostructures and high-temperature cuprate superconductors [41, 42, 43].
The advantage of HubbardNet is that the fermionic wavefunction symmetry can be
directly incorporated into the model hamiltonian and the basis set
construction, similar to ED. For example, in twisted bilayer transition metal
dichalcogenides (TMDCs), where unconventional superconductivity has been
reported [44], the Coulomb interaction $U$ strength is related to the band
flatness and $t$ is related to the bandwidth, both of which are functions of
the twist angle $\theta$. It has been shown both theoretically [4, 45] and
experimentally [46] that twisted bilayer TMD can be accurately represented
with a Fermi-Hubbard model on a triangular, honeycomb, or Kagome lattice. From
the results of the first-principles calculations, one could derive reduced-
band models to obtain the $U/t$ ratio as a function of $\theta$. HubbardNet
essentially allows for the mapping of the phase boundary in twisted TMDCs as a
function of $\theta$ with a single training. Our code is openly available
[47].

###### Acknowledgements.

Z.Z. and E.K. acknowledge the support from the STC Center for Integrated
Quantum Materials, NSF Grant No. DMR-1231319, and NSF DMREF Grant No. 1922172,
and a Simons Foundation Award No. 896626. Calculations were performed on the
Cannon cluster supported by the FAS Division of Science, Research Computing
Group at Harvard University, Google Colaboratory, and the Sherlock Cluster
supported by the Stanford Research Computing Center.

## References

  * Anderson [1998] P. Anderson, The theory of superconductivity in the high t, Tc (1998).
  * Huang  _et al._ [2017] E. W. Huang, C. B. Mendl, S. Liu, S. Johnston, H.-C. Jiang, B. Moritz, and T. P. Devereaux, Numerical evidence of fluctuating stripes in the normal state of high-Tc cuprate superconductors, Science 358, 1161 (2017), arXiv:1612.05211 [cond-mat.str-el] .
  * Bedell  _et al._ [1990] K. S. Bedell, D. Coffey, D. Meltzer, D. Pines, and J. Schrieffer,  _High temperature superconductivity: Proceedings_ , Tech. Rep. (Los Alamos National Lab., NM (USA), 1990).
  * Wu  _et al._ [2018] F. Wu, T. Lovorn, E. Tutuc, and A. H. MacDonald, Hubbard model physics in transition metal dichalcogenide moiré bands, Phys. Rev. Lett. 121, 026402 (2018).
  * Xie  _et al._ [2019] Y. Xie, B. Lian, B. Jäck, X. Liu, C.-L. Chiu, K. Watanabe, T. Taniguchi, B. A. Bernevig, and A. Yazdani, Spectroscopic signatures of many-body correlations in magic-angle twisted bilayer graphene, Nature (London) 572, 101 (2019), arXiv:1906.09274 [cond-mat.mes-hall] .
  * Pan  _et al._ [2020a] H. Pan, F. Wu, and S. Das Sarma, Band topology, hubbard model, heisenberg model, and dzyaloshinskii-moriya interaction in twisted bilayer ${\mathrm{wse}}_{2}$, Phys. Rev. Research 2, 033087 (2020a).
  * Pan  _et al._ [2020b] H. Pan, F. Wu, and S. Das Sarma, Quantum phase diagram of a moiré-hubbard model, Phys. Rev. B 102, 201104 (2020b).
  * Ochi  _et al._ [2018] M. Ochi, M. Koshino, and K. Kuroki, Possible correlated insulating states in magic-angle twisted bilayer graphene under strongly competing interactions, Phys. Rev. B 98, 081102 (2018).
  * Dodaro  _et al._ [2018] J. F. Dodaro, S. A. Kivelson, Y. Schattner, X. Q. Sun, and C. Wang, Phases of a phenomenological model of twisted bilayer graphene, Phys. Rev. B 98, 075154 (2018), arXiv:1804.03162 [cond-mat.supr-con] .
  * Wong  _et al._ [2020] D. Wong, K. P. Nuckolls, M. Oh, B. Lian, Y. Xie, S. Jeon, K. Watanabe, T. Taniguchi, B. A. Bernevig, and A. Yazdani, Cascade of electronic transitions in magic-angle twisted bilayer graphene, Nature (London) 582, 198 (2020), arXiv:1912.06145 [cond-mat.mes-hall] .
  * Potasz  _et al._ [2021] P. Potasz, M. Xie, and A. H. MacDonald, Exact Diagonalization for Magic-Angle Twisted Bilayer Graphene, Phys. Rev. Lett. 127, 147203 (2021), arXiv:2102.02256 [cond-mat.str-el] .
  * Morales-Durán  _et al._ [2021] N. Morales-Durán, A. H. MacDonald, and P. Potasz, Metal-insulator transition in transition metal dichalcogenide heterobilayer moiré superlattices, Phys. Rev. B 103, L241110 (2021), arXiv:2011.13558 [cond-mat.str-el] .
  * Scalapino  _et al._ [1992] D. J. Scalapino, S. R. White, and S. C. Zhang, Superfluid density and the drude weight of the hubbard model, Phys. Rev. Lett. 68, 2830 (1992).
  * Freericks and Monien [1994] J. K. Freericks and H. Monien, Phase diagram of the Bose-Hubbard model, Europhysics Letters (EPL) 26, 545 (1994).
  * Kühner and Monien [1998] T. D. Kühner and H. Monien, Phases of the one-dimensional Bose-Hubbard model, Phys. Rev. B 58, R14741 (1998).
  * Capogrosso-Sansone  _et al._ [2007] B. Capogrosso-Sansone, N. V. Prokof’ev, and B. V. Svistunov, Phase diagram and thermodynamics of the three-dimensional Bose-Hubbard model, Phys. Rev. B 75, 134302 (2007).
  * Carleo and Troyer [2017] G. Carleo and M. Troyer, Solving the quantum many-body problem with artificial neural networks, Science 355, 602 (2017), arXiv:1606.02318 [cond-mat.dis-nn] .
  * Nomura  _et al._ [2017] Y. Nomura, A. S. Darmawan, Y. Yamaji, and M. Imada, Restricted boltzmann machine learning for solving strongly correlated quantum systems, Phys. Rev. B 96, 205152 (2017).
  * Bravyi  _et al._ [2019] S. Bravyi, D. Gosset, R. König, and K. Temme, Approximation algorithms for quantum many-body problems, Journal of Mathematical Physics 60, 032203 (2019), arXiv:1808.01734 [quant-ph] .
  * Saito [2017] H. Saito, Solving the Bose-Hubbard Model with Machine Learning, Journal of the Physical Society of Japan 86, 093001 (2017), arXiv:1707.09723 [cond-mat.dis-nn] .
  * Saito and Kato [2018] H. Saito and M. Kato, Machine learning technique to find quantum many-body ground states of bosons on a lattice, Journal of the Physical Society of Japan 87, 014001 (2018), arXiv:1709.05468 [cond-mat.dis-nn] .
  * McBrian  _et al._ [2019] K. McBrian, G. Carleo, and E. Khatami, Ground state phase diagram of the one-dimensional Bose-Hubbard model from restricted Boltzmann machines, in  _Journal of Physics Conference Series_, Journal of Physics Conference Series, Vol. 1290 (2019) p. 012005, arXiv:1903.03076 [cond-mat.str-el] .
  * Choo  _et al._ [2020] K. Choo, A. Mezzacapo, and G. Carleo, Fermionic neural-network states for ab-initio electronic structure, Nature Communications 11, 2368 (2020), arXiv:1909.12852 [physics.comp-ph] .
  * Robledo Moreno  _et al._ [2021] J. Robledo Moreno, G. Carleo, A. Georges, and J. Stokes, Fermionic Wave Functions from Neural-Network Constrained Hidden States, arXiv e-prints , arXiv:2111.10420 (2021), arXiv:2111.10420 [cond-mat.str-el] .
  * Yoshioka  _et al._ [2021] N. Yoshioka, W. Mizukami, and F. Nori, Solving quasiparticle band spectra of real solids using neural-network quantum states, Communications Physics 4, 106 (2021), arXiv:2010.01358 [cond-mat.str-el] .
  * Bennewitz  _et al._ [2021] E. R. Bennewitz, F. Hopfmueller, B. Kulchytskyy, J. Carrasquilla, and P. Ronagh, Neural Error Mitigation of Near-Term Quantum Simulations, arXiv e-prints , arXiv:2105.08086 (2021), arXiv:2105.08086 [quant-ph] .
  * Hartmann and Carleo [2019] M. J. Hartmann and G. Carleo, Neural-Network Approach to Dissipative Quantum Many-Body Dynamics, Phys. Rev. Lett. 122, 250502 (2019), arXiv:1902.05131 [quant-ph] .
  * Nagy and Savona [2019] A. Nagy and V. Savona, Variational Quantum Monte Carlo Method with a Neural-Network Ansatz for Open Quantum Systems, Phys. Rev. Lett. 122, 250501 (2019), arXiv:1902.09483 [quant-ph] .
  * Vicentini  _et al._ [2019] F. Vicentini, A. Biella, N. Regnault, and C. Ciuti, Variational Neural-Network Ansatz for Steady States in Open Quantum Systems, Phys. Rev. Lett. 122, 250503 (2019), arXiv:1902.10104 [quant-ph] .
  * Ponomarev  _et al._ [2009] A. V. Ponomarev, S. Denisov, and P. Hänggi, ac-Driven Atomic Quantum Motor, Phys. Rev. Lett. 102, 230601 (2009), arXiv:0902.0489 [physics.atom-ph] .
  * Ponomarev  _et al._ [2010] A. V. Ponomarev, S. Denisov, and P. Hänggi, Lévy distribution in many-particle quantum systems, Phys. Rev. A 81, 043615 (2010), arXiv:0907.4328 [cond-mat.quant-gas] .
  * Raventós  _et al._ [2017] D. Raventós, T. Graß, M. Lewenstein, and B. Juliá-Díaz, Cold bosons in optical lattices: a tutorial for exact diagonalization, Journal of Physics B Atomic Molecular Physics 50, 113001 (2017), arXiv:1410.7280 [cond-mat.quant-gas] .
  * Blankenbecler  _et al._ [1981] R. Blankenbecler, D. J. Scalapino, and R. L. Sugar, Monte carlo calculations of coupled boson-fermion systems. i, Phys. Rev. D 24, 2278 (1981).
  * Loshchilov and Hutter [2016] I. Loshchilov and F. Hutter, SGDR: Stochastic Gradient Descent with Warm Restarts, arXiv e-prints , arXiv:1608.03983 (2016), arXiv:1608.03983 [cs.LG] .
  * Fisher  _et al._ [1989] M. P. A. Fisher, P. B. Weichman, G. Grinstein, and D. S. Fisher, Boson localization and the superfluid-insulator transition, Phys. Rev. B 40, 546 (1989).
  * Bloch [2005] I. Bloch, Ultracold quantum gases in optical lattices, Nature Physics 1, 23 (2005).
  * Ejima  _et al._ [2012] S. Ejima, H. Fehske, F. Gebhard, K. zu Münster, M. Knap, E. Arrigoni, and W. von der Linden, Characterization of mott-insulating and superfluid phases in the one-dimensional Bose-Hubbard model, Phys. Rev. A 85, 053644 (2012).
  * Kashurnikov  _et al._ [1996] V. A. Kashurnikov, A. V. Krasavin, and B. V. Svistunov, Mott-insulator-superfluid-liquid transition in a one-dimensional bosonic Hubbard model: Quantum Monte Carlo method, Soviet Journal of Experimental and Theoretical Physics Letters 64, 99 (1996).
  * Bose Hubbard  _et al._ [2011] S. Bose Hubbard, H. Fehske, and F. Gebhard, Dynamic properties of the one-dimensional Bose-Hubbard model, Europhysics Letters 93, 30002 (2011).
  * Carleo  _et al._ [2019] G. Carleo, K. Choo, D. Hofmann, J. E. T. Smith, T. Westerhout, F. Alet, E. J. Davis, S. Efthymiou, I. Glasser, S.-H. Lin, M. Mauri, G. Mazzola, C. B. Mendl, E. van Nieuwenburg, O. O’Reilly, H. Théveniaut, G. Torlai, F. Vicentini, and A. Wietek, Netket: A machine learning toolkit for many-body quantum systems, SoftwareX , 100311 (2019).
  * Cao  _et al._ [2018] Y. Cao, V. Fatemi, S. Fang, K. Watanabe, T. Taniguchi, E. Kaxiras, and P. Jarillo-Herrero, Unconventional superconductivity in magic-angle graphene superlattices, Nature (London) 556, 43 (2018), arXiv:1803.02342 [cond-mat.mes-hall] .
  * Fradkin  _et al._ [2015] E. Fradkin, S. A. Kivelson, and J. M. Tranquada, Colloquium: Theory of intertwined orders in high temperature superconductors, Reviews of Modern Physics 87, 457 (2015), arXiv:1407.4480 [cond-mat.supr-con] .
  * Heidrich-Meisner  _et al._ [2004] F. Heidrich-Meisner, A. Honecker, D. C. Cabra, and W. Brenig, Comment on “anomalous thermal conductivity of frustrated heisenberg spin chains and ladders”, Phys. Rev. Lett. 92, 069703 (2004).
  * Wang  _et al._ [2020] L. Wang, E.-M. Shih, A. Ghiotto, L. Xian, D. A. Rhodes, C. Tan, M. Claassen, D. M. Kennes, Y. Bai, B. Kim, K. Watanabe, T. Taniguchi, X. Zhu, J. Hone, A. Rubio, A. N. Pasupathy, and C. R. Dean, Correlated electronic phases in twisted bilayer transition metal dichalcogenides, Nature Materials 19, 861 (2020).
  * Angeli  _et al._ [2022] M. Angeli, G. R. Schleder, and E. Kaxiras, Twistronics of Janus transition metal dichalcogenide bilayers, arXiv e-prints , arXiv:2207.05788 (2022), arXiv:2207.05788 [cond-mat.mtrl-sci] .
  * Xu  _et al._ [2022] Y. Xu, K. Kang, K. Watanabe, T. Taniguchi, K. F. Mak, and J. Shan, A tunable bilayer Hubbard model in twisted WSe2, Nature Nanotechnology 17, 934 (2022), arXiv:2202.02055 [cond-mat.str-el] .
  * Zhu  _et al._ [2022] Z. Zhu, M. Mattheakis, W. Pan, and E. Kaxiras, Model for HubbardNet: Efficient Predictions of the Bose-Hubbard Model Spectrum with Deep Neural Networks: https://github.com/ziyanzzhu/HubbardNet (2022).


00footnotetext: The second author partially supported by the Israel Science
Foundation,grantno.1970/18

# Major index on involutions

Eli Bagno  Department of Applied Mathematics, Jerusalem College of Technology,
21 Havaad Haleumi St., Jerusalem, Israel<EMAIL_ADDRESS>and  Yisca Kares
Department of Mathematics, Bar-Ilan University, Ramat-Gan 52900, Israel
<EMAIL_ADDRESS>
###### Abstract.

We find the range of the major index on the various conjugacy classes of
involutions in the symmetric group $S_{n}$. In addition to indicating the
minimum and the maximum values, we show that except for the case of
involutions without fixed points, all the values in the range are attained.
For the conjugacy classes of involutions without fixed points, we show that
the only missing values are one more than the minimum and one less than the
maximum.

##  1\. Introduction

Let ${\mathcal{S}}_{n}$ be the symmetric group on $n$ elements, let
$\lambda=(\lambda_{1},\dots,\lambda_{k})$ be a partition of $n$ which will be
identified with its Young diagram of the shape $\lambda$, and let
$SYT(\lambda)$ denote the set of all standard Young tableaux of shape
$\lambda$. We say that $i\in\\{1,\dots,n-1\\}$ is a descent of a permutation
$\pi\in{\mathcal{S}}_{n}$ if $\pi(i)>\pi(i+1)$.

The major index of a permutation $\pi\in{\mathcal{S}}_{n}$ , $maj(\pi)$, is
the sum of the descents of $T$. Likewise, if $T\in SYT(\lambda)$ is a standard
tableaux of shape $\lambda$ then $i\in\\{1,\dots,n-1\\}$ is a descent of $T$
if $i+1$ is located in $T$ in a row lower than $i$. Again, the major index,
$maj(T)$, is the sum of the descents. See a more elaborated exposition of
these concepts in Section 2.
Owing to the RSK algorithm which associates a pair of standard tableaux of the
same shape to every permutation $\pi\in{\mathcal{S}}_{n}$, these two versions
of the major index are tightly connected.

The generating function of the major index over the entire group
${\mathcal{S}}_{n}$ is known since the days of McMahon [3] who showed that the
distribution of the major index on all permutations in ${\mathcal{S}}_{n}$ is
the same as the distribution of inversions. Explicitly:

$\sum\limits_{\pi\in{\mathcal{S}}_{n}}q^{maj(\pi)}=[n]_{q}!$

where $[n]_{q}=\frac{q^{n}-1}{q-1}$ and
$[n]_{q}!=[n]_{q}\cdot[n-1]_{q}\cdots[1]_{q}$.

On the other hand, the generating function

$SYT(\lambda)^{maj}(q):=\sum\limits_{T\in SYT(\lambda)}q^{maj(T)}=\sum_{k\geq
0}m_{\lambda,k}q^{k}$

has two elegant closed forms: one due to Steinberg [8], based on dimensions of
irreducible representations of $GL_{n}(\mathbb{F}_{q})$, and the other is
based on Stanely’s generalization of the hook formula [6]. It reads:

$SYT(\lambda)^{maj}(q)=\frac{q^{b(\lambda)}[n]_{q}!}{\Pi_{c\in\lambda}[h_{c}]_{q}}$

where $h_{c}$ is the hook length of the cell $c$ and
$b(\lambda)=\sum\limits_{i>0}(i-1)\lambda_{i}$. The sequence
$\\{m_{\lambda,k}\\}$ is called the fake degree sequence and has appeared in
various algebraic and representation-theoretic contexts, such as the degree
polynomials of unipotent $GL_{n}(\mathbb{F}_{q})$ -representations due to
Green [2]. Another context in which the fake degrees have a significant rule
is the coinvariant algebra of ${\mathcal{S}}_{n}$. More broadly, the $k$-th
part of the coinvariant algebra of ${\mathcal{S}}_{n}$ decomposes into
irreducible representations, each of them appearing with multiplicity
$m_{k,\lambda}$ .

The range of the major index inside a single shape has been recently explored
by Billey, Konvalinka and Swanson in [1] (see Theorem 1.1. there). We
concentrate here on the range of values the major index attains at the
conjugacy classes of involutions (i.e. permutations of order less or equal to
$2$). This will be done using the RSK correspondence between permutations and
pairs of standard Young tableaux, which turns out to be a bijection with
standard Young tableaux when restricted to involutions.

Recall that conjugacy classes of involutions are determined by the number of
fixed points. According to a result by Schützenberger [7], all the involutions
of a given conjugacy class with $r$ fixed points are associated by RSK with
standard tableaux of shapes having $r$ odd columns.

The combination of these two results, together with an algorithm that runs
through the standard Young tableaux of a fixed number of odd columns enable us
to calculate the range of the major index over the entire conjugacy classes.

The paper is structured as follows: First, we provide some background in the
following section, including the definitions of the statistics we work with,
the Young diagrams, the Young tableaux, and the Robinson-Schensted
correspondence. Additionally, we present some preparatory lemmas in Section 2
to enhance the readability of our proof of the main theorem. Section 3
introduces our main theorem and the algorithm we use to demonstrate it, along
with the proof’s overarching idea. Further, we justify the algorithm’s steps
in Section 4, where we provide detailed proofs for each step. Finally, we
cover some special cases at the end of Section 4.

##  2\. Background

###  2.1. Statistics on permutations

We define here two of the well-known statistics on ${\mathcal{S}}_{n}$,
namely, the descent number and the major index.

######  Definition 2.1.

Let $\pi\in{\mathcal{S}}_{n}$. The descent set of $\pi$ is:

$Des(\pi)=\\{i\in\\{1,\dots,n-1\\}\mid\pi(i)>\pi(i+1)\\}.$

We also define $des(\pi)=|Des(\pi)|$.

The major index of $\pi$ is defined by:

$maj(\pi)=\sum\limits_{i\in Des(\pi)}{i}.$

######  Example 2.2.

Let $\pi=5321476\in\mathcal{S}_{7}$. Then $Des(\pi)=\\{1,2,3,6\\}$,
$des(\pi)=4$ and $maj(\pi)=1+2+3+6=12$.

###  2.2. Young diagrams and standard tableaux

A Young diagram is a finite collection of cells in the plane, arranged in
left-justified rows, such that the row lengths are increasing. The sequence
listing the numbers of cells in each row gives a partition $\lambda$ of a non-
negative integer $n$. The Young diagram is said to be of shape $\lambda$.
Also, let $\lambda^{\prime}$ be the transpose of $\lambda$. We often use Greek
letters to denote the diagrams as well as their shapes (See figure 1 for an
example).

$\lambda=$ &

$\lambda^{\prime}=$ &

Figure 1.  $\lambda=(4,3,1,1),\lambda^{\prime}=(4,2,2,1)$

A Young tableau is a filling of a young diagram a by the numbers $1,\dots,n$.
A tableau is standard if the numbers are increasing through rows and column.
(see example in figure 2).

1 & 4 5 6

2 7 9

3

8

Figure 2.  A standard tableau $T$ of shape $\lambda=(4,3,1,1)$,
$Des(T)=\\{1,2,6,7\\},maj(T)=1+2+6+7=16$

For a shape $\lambda$, let $\operatorname{SYT}(\lambda)$ be the set of
standard Young tableaux of shape $\lambda$. The descent set of $T$ is

$\operatorname{Des}(T):=\\{i\,\mid i+1\text{ appears in a lower row of
}T\text{ than }i\\}.$

Define also the major index of a standard Young tableau $T$ by

$maj(T)=\sum_{i\in Des(T)}{i}.$

(See example in figure 2).

###  2.3. The RSK algorithm

We shall make use of the Robinson-Schensted-Knuth (RSK) correspondence which
maps each permutation $\pi\in{\mathcal{S}}_{n}$ to a pair $(P_{\pi},Q_{\pi})$
of standard Young tableaux of the same shape $\lambda$. A detailed description
can be found, for example, in [4, Ch. 3.1] or in [5, Ch. 7.11].

The following is a crucial property of the RSK correspondence which will be
vastly used in this paper.

######  Fact 2.3.

For each $\pi\in{\mathcal{S}}_{n}$ one has

$Q_{\pi}=P_{\pi^{-1}}.$

Note also that the RSK correspondence is a $\operatorname{Des}$-preserving and
hence also a $maj$\- preserving bijection in the following sense:

######  Fact 2.4.

For every permutation $\pi\in{\mathcal{S}}_{n}$,

$\operatorname{Des}(P_{\pi})=\operatorname{Des}(\pi^{-1})\quad\text{ and
}\quad\operatorname{Des}(Q_{\pi})=\operatorname{Des}(\pi).$

Consequently, we have also $maj(P_{\pi})=maj(\pi^{-1})$ as well as
$maj(Q_{\pi})=maj(\pi)$.

It follows from Fact 2.3 that $\pi$ is an involution if and only if
$P_{\pi}=Q_{\pi}$ so that by restricting the RSK correspondence to the set of
involutions which will be denoted here by $I_{n}$, we get a $Des-$ preserving
bijection from $I_{n}$ to the set of standard Young tableaux of order $n$,
$SYT(n)$.

###  2.4. Conjugacy classes of involutions

Conjugacy classes in ${\mathcal{S}}_{n}$ are determined by their cycle
structures, which are partitions of $n$. The conjugacy classes of involutions
in $I_{n}$ are of cycle structure $(2^{k},1^{r})$ (where $a^{b}$ means $b$
appearances of $a$) such that $0\leq r\leq n$, $0\leq k\leq\frac{n}{2}$ and
$2k+r=n$. In other words, conjugacy classes of involutions are determined by
the number of fixed points.

The following known result by Schützenberger [7] gives a full description of
the image of each conjugacy class of involutions under the RSK correspondence.

######  Proposition 2.5.

An involution $\pi\in I_{n}$ has $r$ fixed points if and only if $P(\pi)$ has
$r$ columns of odd length.

In light of this characterization, we present the following notations:

######  Definition 2.6.

The set of Young diagrams of size $n$ having exactly $r$ odd columns will be
denoted by $D_{n}(r)$. The set of standard Young tableaux of shapes taken from
$D_{n}(r)$ will be denoted by $SYT_{n}(r)$.

The discussion above can now be concisely formulated as follows:

######  Proposition 2.7.

Let ${\mathcal{C}}_{\mu}$ be the conjugacy class corresponding to the
partition $\mu=(2^{k},1^{r})$. Then the restriction of the RSK correspondence

$R:{\mathcal{C}}_{\mu}\rightarrow SYT_{n}(r)$

is a bijection that preserves the major index, i.e. for each
$\pi\in{\mathcal{C}}_{\mu}$ we have $maj(\pi)=maj(R(\pi))$.

The following definition will be needed in the sequel (note the slight change
from a similar definition made in the introduction; we prefer this notation
for the sake of convenience).

######  Definition 2.8.

For a shape $\lambda=(\lambda_{0},\lambda_{1},\dots,\lambda_{u})$, let

$b(\lambda)=\sum\limits_{i=0}^{u}{i\lambda_{i}}.$

######  Remark 2.9.

The numbers $b(\lambda)$, ($b(\lambda^{\prime})$) can be easily calculated by
writing for each $i$ the number $i$ inside each square of row (column) $i$ of
$\lambda$ and adding up the numbers to get the rows (columns) sums
respectively.

0 &0 0 0
1 1 1
2 2 , 0 &1 2 3
0 1 2
0 1

Figure 3. Calculating $b(\lambda)$ (left) and $b(\lambda^{\prime})$ (right) by
adding up rows and columns respectively

.

The continuity of $maj$ inside the Young tableaux of a single Young diagram of
a specific shape has been recently proven by Billey, Konvalinka, and Swanson
in [1] (see Theorem 1.1. there). The following is a reformulation of their
result:

######  Proposition 2.10.

Let $\lambda$ be a Young diagram. Then we have:

$m(\lambda):={\rm Min}\\{maj(T)\mid T\in SYT(\lambda)\\}=b(\lambda).$

$M(\lambda):={\rm Max}\\{maj(T)\mid T\in SYT(\lambda)\\}={n\choose
2}-b(\lambda^{\prime}).$

Moreover, every value between $m(\lambda)$ and $M(\lambda)$ appears at least
once, except when $\lambda$ is a rectangle with at least two rows and columns,
in which case the values $m(\lambda)+1$ and $M(\lambda)-1$ are missing.

The following lemma determines the diagrams of $D_{n}(r)$ which attain the
minimum and the maximum of the major index.

######  Lemma 2.11.

Let $n=2k+r$.

  1. (1)

The minimum value of the major index on $D_{n}(r)$ is $k$. It is attained by
the diagram $\lambda=(n-k,k)$.

  2. (2)

The maximum value of the major index on $D_{n}(r)$ is ${n\choose 2}-{r\choose
2}$. It is attained by the diagram $\lambda=(r,1^{2k})$.

###### Proof.

  1. (1)

Let $m={\rm Min}\\{maj(T)\mid T\in SYT_{n}(r)\\}$ and let $\lambda=(n-k,k)$.
Since $n=2k+r$, it is easy to see that $\lambda$ has $r$ odd columns, so that
$\lambda\in D_{n}(r)$. By Proposition 2.10, $b(\lambda)=0\cdot(n-k)+1\cdot
k=k$ is the minimum value of the major index on $\lambda$, hence $m\leq k$.
Now, for each $\nu\in D_{n}(r)$, if $b(\nu)<k$ then the complement of the
first line of $\nu$ contains less than $k$ squares and so the first line of
$\nu$ must contain more than $n-k=r+2k-k=r+k$ squares. As the number of odd
columns is $r$, there are more than $k$ squares in the first line which start
columns of size at least $2$ so that $b(\nu)>k$, a contradiction.

  2. (2)

Let $M={\rm Max}\\{maj(T)\mid T\in SYT_{n}(r)\\}$ and let
$\lambda=(r,1^{2k})$. It is easy to see that $\lambda\in D_{n}(r)$. We have
$b(\lambda^{\prime})=1+2+\cdots+(r-1)={r\choose 2}$ so we conclude from
Proposition 2.10 that $M\geq{n\choose 2}-{r\choose 2}$. Every $\nu\in
D_{n}(r)$ must have at least $r$ columns, hence we must have
$b(\nu^{\prime})\geq{r\choose 2}$, and we are done.

∎

We will show in the sequel that for most of the conjugacy classes of
involutions, the range of the major index on these classes is an interval, but
we start with the exceptional case: the class of involutions without fixed
points.

######  Lemma 2.12.

Let ${\mathcal{C}}$ be the conjugacy class corresponding to the partition
$\mu=(2^{k})$. Then:

  1. (1)

There is no $\pi\in{\mathcal{C}}$ such that $maj(\pi)=k+1$.

  2. (2)

There is no $\pi\in{\mathcal{C}}$ such that $maj(\pi)={n\choose 2}-1$.

###### Proof.

According to Proposition 2.7, it is sufficient to prove that there is no $T\in
SYT_{n}(0)$ such that $maj(T)\in\\{k+1,{n\choose 2}-1\\}$.

  1. (1)

First, since there are no odd columns in any diagram of $D_{n}(0)$ and $n=2k$,
the number of columns does not exceed $k$. If there is some $T\in SYT_{n}(0)$
such that $maj(T)=k+1$, then we have either $Des(T)=\\{k+1\\}$ or $Des(T)$
contains elements less than $k+1$, the sum of which is $k+1$. In the first
case, all the elements $1,\dots,k+1$ occupy the first line, which means that
$T$ has more than $k$ columns, a contradiction. In the second case,
$k+1,\dots,2k$ must be placed in different columns. Since all the columns
contain at least $2$ elements, the numbers $1,\dots,k$ are also placed in
different columns but this means that the major index is $k$ and not $k+1$.

  2. (2)

We have already seen in Lemma 2.11 that the maximal value, $M={n\choose 2}$ is
attained by the one-column diagram. Any other diagram in $D_{n}(0)$ contains
at least two columns, all of them even. If $T\in SYT_{n}(0)$ is such that
$maj(T)={n\choose 2}-1=2+\cdots+(n-1)$, then we must have
$Des(T)=\\{2,\dots,n-1\\}$. The element $1$ must be placed in the upper left
corner and since $1\notin Des(T)$, the element $2$ must be placed in the first
line, right after the element $1$. Since $2\in Des(T)$, we must put $3$ below
$1$ and this argument proceeds further and gives us the tableau of Fig. 4
which in not in $D_{n}(0)$.

$\lambda=$ 1 & 2

3

4

⋅

⋅

⋅

n

Figure 4.

∎

######  Definition 2.13.

Recall that a diagram of the form $(u,1^{n-u})$ is called a hook. If the
length of the leg of $\lambda$, $n-u+1$ is odd then $\lambda$ will be called
an odd hook. It will be called an even hook otherwise.

The following lemma is a direct consequence of Proposition 2.7 and Lemma
2.11.2.

######  Lemma 2.14.

Let $\mu$ be the conjugacy class of involutions with exactly $r$ fixed points.
Then the maximal value of major index over $\mu$ is the maximal value of
$D_{n}(r)$ which is obtained by the hook $\lambda=(r,1^{n-r})$.

The following lemma will be of great use when we prove our main result in the
next section.

######  Lemma 2.15.

Let $n\geq 6$ and let $\lambda\vdash n$ be such that $\lambda\neq(1^{n})$ and
$\lambda\neq(n)$. Then $M(\lambda)-m(\lambda)\geq 4$.

Moreover, if $\lambda=(a^{b})$ is a rectangle then $M(\lambda)-m(\lambda)\geq
6$.

###### Proof.

Recall from lemma 2.10 that $m(\lambda)=b(\lambda)$, while
$M(\lambda)={n\choose 2}-b(\lambda^{\prime})$, so that
$\Delta:=M(\lambda)-m(\lambda)={n\choose 2}-(b(\lambda)+b(\lambda^{\prime}))$.

Let $T_{1}$ be the tableau of shape $\lambda$, filled by the numbers
$0,\dots,n-1$ in such a way that they can be read in this order along the
consecutive rows, increasing along each row.

Let $T_{2}$ be the tableau of shape $\lambda$, filled in such a way that the
square placed in row $i$ and column $j$ contains $i+j$ where the rows and
columns are numbered $0,1,\dots$, etc.

Here is an example with $\lambda=(4,3,1,1)$:

$T_{1}=$ boxsize=3em 0 & 1 2 3
4 5 6
7
8
boxsize=3em $T_{2}=$ 0+0 & 0+1 0+2 0+3
1+0 1+1 1+2
2+0
3+0
boxsize=normal

Note that the sum of elements of $T_{1}$ is ${n\choose 2}$ while the sum of
elements of $T_{2}$ is $b(\lambda)+b(\lambda^{\prime})$, so that $\Delta$ is
the difference of the row sums. Denote for each $i$ by $\Delta_{i}$ the
difference between the sum of elements of row number $i$ in $T_{1}$ and the
sum of elements of row number $i$ in $T_{2}$, so that
$\Delta=\sum\limits_{i}\Delta_{i}$.

It is easy to see that $\Delta_{0}=0$ and
$\Delta_{1}=\lambda_{1}(\lambda_{0}-1)$. Moreover, for each $i\geq 2$, we have
$\Delta_{i}\geq 0$ since the first element of each such row in $T_{1}$ is
greater than the corresponding element in $T_{2}$ and in both tableaux all
elements of row $i$ are consecutive.

By the assumptions of the lemma, we must have $\lambda_{1}\geq 1$. we divide
in several cases according to the value of $\lambda_{1}$:

  * •

If $\lambda_{1}\geq 3$ then $\lambda_{0}-1\geq 2$ so that
$\Delta\geq\Delta_{1}\geq 6$ and we are done.

  * •

If $\lambda_{1}=2$ then we must have $\lambda_{0}\geq 2$. Now, if
$\lambda_{0}=2$ then $\lambda$ contains the diagram $(2,2,1)$ which satisfies
$\Delta\geq 4$ by inspection.

The case $\lambda_{1}=2$ and $\lambda_{0}=3$ is symmetric to the last case,
i.e $\Delta$ is invariant under transposition. Furthermore, if
$\lambda_{0}\geq 4$ then $\Delta_{1}\geq 6$ and we are done.

  * •

if $\lambda_{1}=1$ then we must have $\lambda_{0}\geq 2$ since otherwise we
get $\lambda=(1^{n})$ which contradicts the assumptions. If $\lambda_{0}=2$
then $\lambda$ contains the diagram $(2,1^{4})$, since by the assumption
$n\geq 6$, so that $\Delta\geq 4$ by inspection.

Regarding the second claim, if $\lambda=(a^{b})$ is a rectangle, then since
$n\geq 6$, it must contain either the shape $(3^{2})$ or the shape $(2^{3})$.
If $\lambda$ contains the shape $(3^{2})$ then it is easy to see that
$\Delta\geq\lambda_{1}(\lambda_{0}-1)\geq 6$. The other case follows by
symmetry. ∎

##  3\. The range of the major index on conjugacy classes of involutions

Our main result is the following:

######  Theorem 3.1.

Let $\mu=(2^{k},1^{r})$ be a partition of $n$ and let ${\mathcal{C}}_{\mu}$ be
the corresponding conjugacy class of involutions in ${\mathcal{S}}_{n}$. Then

  * •

If $r\neq 0$, then the major index on ${\mathcal{C}}_{\mu}$ attains all values
between $k$ and $\binom{n}{2}-\binom{r}{2}$.

  * •

If $r=0$, then the major index on ${\mathcal{C}}_{\mu}$ attains all the values
mentioned in the first clause, excluding $k+1$ and $\binom{n}{2}-1$.

Moreover, any other value outside this range is not attained.

We present here the sketch of the proof and postpone the details to Section 4.

###  3.1. Sketch of the proof

By Propositions 2.5 and 2.7, it is sufficient to prove our results for the
sets $SYT_{r}(n)$, consisting of all standard tableaux of shapes having
exactly $r$ odd columns.

This will be done by ordering the set $D_{n}(r)$ of diagrams of size $n$ with
exactly $r$ odd columns, (see Definition 2.6) by reverse dominant order (see
Definition 2.2.2 in [4]) and presenting an algorithm which starts with the
diagram $\lambda^{0}=(n-k,k)=(k+r,k)$, attaining the minimum value of $maj$
over $D_{n}(r)$ which is $k$ and ends with the odd hook diagram
$\lambda^{e}=(r,1^{2k})$, attaining the maximum value of $maj$ over $D_{n}(r)$
which is ${n\choose 2}-{r\choose 2}$ (see Lemma 2.11). The algorithm traverses
the set $D_{n}(r)$ in such a way that in each step one or two squares of a
diagram $\lambda\in D_{n}(r)$ are transferred to a new place to obtain another
diagram $\nu\in D_{n}(r)$ such that following condition is satisfied:

$M(\lambda)\geq m(\nu),$ (1)

where $M(\lambda)$ ($m(\nu)$) is the maximum (minimum) value of $maj$ on
$SYT(\lambda)$ ($SYT(\nu)$), respectively as in Proposition 2.10. Together
with Proposition 2.10, we are done, apart from some special cases that will be
treated separately in Section 4.

###  3.2. The algorithm

Let $\lambda=(\lambda_{0},\lambda_{1},\dots,\lambda_{t})$. We add infinite
number of zeroes at the end of $\lambda$ and write $\lambda$ as a Young
diagram. Also, for each $i$, denote the last square of the row $\lambda_{i}$
by $\lambda_{i}^{*}$ . Now perform the following steps:

  1. (1)

If $\lambda$ is an odd hook, then we are done by Lemma 2.11.2.

  2. (2)

If $\lambda$ is an even hook, then we distinguish between two cases:

     * •

If $\lambda=(1^{2k})$, i.e. $r=0$, then again we are done by Lemma 2.11.2.

     * •

Otherwise, let $\nu$ be the shape obtained from $\lambda$ by removing the
square $\lambda_{0}^{*}$ and placing it at the end of the first column of
$\lambda$. This $\nu$ is an odd hook so we are back in (1) (see the passage
from $\lambda^{7}$ to $\lambda^{8}$ in Example 3.2). We will justify this step
in Lemma 4.1

  3. (3)

If there is some $0\leq j\leq t$ such that
$\lambda_{j}>\lambda_{j+1}>\lambda_{j+2}$, then let $i$ be maximal with
respect to this property. Since $\lambda_{i}>\lambda_{i+1}>\lambda_{i+2}$, we
have in $\lambda$ a column of length $i+1$ and a column of length $i+2$ (see
Fig. 5). Now, remove $\lambda_{i}^{*}$ and place it as the last square of row
$i+2$ and let $\nu$ be the resulting shape. We will elaborate on the utility
of this step in Lemma 4.2

This case is illustrated in Figs. 5 and 6.

$\lambda=$ & *

$\longrightarrow\nu=$ &

*

Figure 5.  $\lambda=(4,3,1,1)$ and $\nu=(3,3,2,1),i=0$.

$\lambda=$ &

*

$\longrightarrow\nu=$ &

*

Figure 6.  $\lambda=(3,2,1)$ and $\nu=(3,1,1,1),i=1$.

  4. (4)

Otherwise, if there is no $0\leq j\leq t$ such that
$\lambda_{j}>\lambda_{j+1}>\lambda_{j+2}$, then there must exist some $j$ such
that $\lambda_{j}>\lambda_{j+1}=\lambda_{j+2}$ (the existence of such $j$ is
guaranteed since we can always choose $j=t$ to get
$\lambda_{t}>\lambda_{t+1}=0=\lambda_{t+2}$). Let $i$ be maximal with respect
to this property and such that $\lambda_{i}>1$ (recall that the case
$\lambda_{i}=\lambda_{j}=1$ for each $i,j$ was treated earlier in cases (1)
and (2)).

For the $i$ we chose, we must have $\lambda_{i-1}=\lambda_{i}$ since
otherwise, we have already treated this case in (3). This means that the
squares $\lambda_{i-1}^{*}$ and $\lambda_{i}^{*}$ form a vertical domino (see
Figures 7 . We distinguish between two sub-cases.

    1. (a)

If we have $\lambda_{i}-\lambda_{i+1}=1$ then, by the maximality of $i$ we
must have $\lambda_{i}=2$. In this case, we put the domino at the end of the
first column (see Figure 7, Left).

    2. (b)

Now, if $\lambda_{i}-\lambda_{i+1}>1$, then we move that domino to the pair of
squares located right after $\lambda_{i+1}^{*}$ and $\lambda_{i+2}^{*}$ (see
Figure 7, Right).

We justify these steps in Lemma 4.3.

  5. (5)

Back to step $(1)$ with $\lambda=\nu$.

&

*

*

$\longrightarrow$ &

*

* &

*

*

$\longrightarrow$ &

*

*

Figure 7. Left: $i=3$, Right: $i=2$

$\lambda^{0}=$ & *

$\rightarrow$ $\lambda^{1}=$ &

*

$\rightarrow$ $\lambda^{2}=$ & *

$\rightarrow$ $\lambda^{3}=$ &

*

$\rightarrow$ $\lambda^{4}=$ & *

*

$\rightarrow$ $\lambda^{5}=$ &

*

$\rightarrow$ $\lambda^{6}=$ &

*

*

$\rightarrow$ $\lambda^{7}=$ & *

$\rightarrow$ $\lambda^{8}=$

Figure 8.

######  Example 3.2.

Let $n=9$ and $\mu=(2^{4},1)$. Then the RSK algorithm takes the elements of
${\mathcal{C}}_{\mu}$ to $D_{9}(1)$ and by Lemma 2.11.1, the diagram which
attains the minimum value of $maj$ is $\lambda^{0}=(5,4)$. The process of the
algorithm ends in diagram $\lambda^{8}$ which is an odd hook that attains the
maximum value of $maj$. The process is depicted in Fig.8. The data of the
tableaux taking part in the process are listed in the table after that figure.

$i$ | $min(\lambda^{i})$ | $max(\lambda^{i})$ | Step used
---|---|---|---
0 | 4 | 20 | 3
1 | 6 | 24 | 3
2 | 8 | 27 | 3
3 | 10 | 29 | 3
4 | 12 | 30 | 4(b)
5 | 16 | 32 | 3
6 | 18 | 33 | 4(a)
7 | 28 | 35 | 2
8 | 36 | 36 | 1

##  4\. Correctness of the algorithm

In this section, we provide comprehensive proof of the correctness of the
algorithm we presented in the previous section. We divide our treatment into
two parts, according to the number of fixed points $r$ in the conjugacy class
of involutions ${\mathcal{C}}_{\mu}$ corresponding to the partition
$\mu=(2^{k},1^{r})$. We start with the case $r\geq 2$ and defer the case $r<2$
to Subsection 4.1.

The main theme of Theorem 3.1 is showing that in each step of the algorithm
when we pass from $\lambda$ to $\nu$, we have $M(\lambda)\geq m(\nu)$. This
inequality will be proved by one of the following methods. If
$\lambda=(n-k,1^{k})$ is an even hook (this happens in step $(2)$ of the
algorithm), then we prove it directly in Lemma 4.1.

Otherwise, we show that one of the following conditions is satisfied.

  1. (1)

$0\leq m(\nu)-m(\lambda)\leq 4$

  2. (2)

$M(\nu)-M(\lambda)=2.$

These conditions, together with Lemma 2.15 imply that whenever we pass from
$\lambda$ to $\nu$ in the process of the algorithm, we have $M(\lambda)\geq
m(\nu)$.

(Indeed, if $0\leq m(\nu)-m(\lambda)\leq 4$ then $m(\nu)\leq m(\lambda)+4$ and
by Lemma 2.15 we have $M(\lambda)-m(\lambda)\geq 4$, so $m(\lambda)+4\leq
M(\lambda)$, hence $m(\nu)\leq M(\lambda)$. The other case is similar, using
Lemma 2.15 for $\nu$ rather than $\lambda$). For more clarification, see Fig.
9.

Figure 9.

Now, in the case where both $\lambda$ and $\nu$ are not rectangles with more
than $2$ rows, the fact that $M(\lambda)\geq m(\nu)$, together with
Proposition 2.10 which assures us that $maj$ obtains every value in their
ranges, guarantee that in the passage from $\lambda$ to $\nu$, all the values
of $maj$ are attained.

If, $\lambda$ is a rectangle with more than $2$ rows, then by Lemma 2.15, we
have that $M(\lambda)-m(\lambda)\geq 6$ and by Proposition 2.10, each value in
the open integral interval $(m(\lambda)+1,M(\lambda)-1)$ is attained, so that
by an argument similar to the one enclosed by the parentheses above, we
immediately get that $m(\nu)\leq M(\lambda)-2$ and we proceed as before. The
case where $\nu$ is a rectangle with more than $2$ rows is treated similarly.

Note that the Algorithm does not necessarily arrive at each shape of
$D_{n}(r)$, but it does prove that all the values from $m(\lambda^{0})$ to
$M(\lambda^{e})$ are obtained by $maj$.

We turn now to explain the validity of steps $(2),(3)$, and $(4)$ of the
algorithm.

######  Lemma 4.1.

If $\lambda=(n-2k+1,1^{2k-1})$ is an even hook and $\nu=(n-2k,1^{2k})$ is the
odd hook we arrive at in step $(2)$ of the algorithm, then $M(\lambda)\geq
m(\nu)$. In other words, step $(2)$ is justified.

###### Proof.

We first note that the number of odd columns has not been changed by the
transfer of the square $\lambda_{0}^{*}$.

Let us write $n=2k+r$, so that $\lambda=(r+1,1^{n-r-1})$ and
$\nu=(r,1^{n-r})$. We have then $M(\lambda)=\binom{n}{2}-\binom{r+1}{2}$ and
$m(\nu)=\binom{n-r+1}{2}$ and we have to show that
$\binom{n}{2}-\binom{r+1}{2}\geq\binom{n-r+1}{2}$.

Showing this inequality is equivalent to showing that $(n-r)^{2}+r^{2}\leq
n^{2}-2n$ which is equivalent in turn to $(n-r-1)(r-1)\geq 1$ .

The last inequality holds if and only if $r\geq 2$, which is the working
assumption of this section ($n-r\geq 2$ is obvious since $\lambda$ is an even
hook). The case $r<2$ will be discussed in Section 4.1.

∎

######  Lemma 4.2.

If we move a square in the diagram $\lambda$ to get the diagram $\nu$
according to the conditions of step $(3)$ of the algorithm, then
$M(\lambda)\geq m(\nu)$, i.e. step (3) is justified.

###### Proof.

The length $i+1$ of the column ending with the square $\lambda^{*}_{i}$ is
reduced by 1, while the length of the column in $\nu$ to which we added the
square, increases from $i+1$ to $i+2$, so that an odd column becomes an even
column and vice versa. There is no change in any other columns.

Note that for each $j\notin\\{i,i+2\\}$ we have $\nu_{j}=\lambda_{j}$, and
also $\nu_{i}=\lambda_{i}-1$ and $\nu_{i+2}=\lambda_{i+2}+1$, so we have

$b(\nu)=i(\lambda_{i}-1)+(i+2)(\lambda_{i+2}+1)+\sum_{j\neq
i,i+2}{j\lambda_{j}}=$ $i\lambda_{i}+(i+2)\lambda_{i+2}+2+\sum_{j\neq
i,i+2}j\lambda_{j}=b(\lambda)+2,$

so that $m(\nu)-m(\lambda)=2\leq 4$. Hence, by the discussion in the beginning
of this section, we conclude that $M(\lambda)\geq m(\nu)$.

∎

######  Lemma 4.3.

If we move a domino in the diagram $\lambda$ to get the diagram $\nu$
according to the conditions of either one of the two cases of step $(4)$ of
the algorithm, then we have $M(\lambda)\geq m(\nu)$, i.e. step $(4)$ is
justified.

###### Proof.

Since we move two squares from one column to another, the number of odd
columns is not changed in the passage from $\lambda$ to $\nu$.

We show now that in case (a) we have $0\leq M(\nu)-M(\lambda)\leq 4$, while in
case (b) we have $0\leq m(\nu)-m(\lambda)\leq 4$ .

  * •

In case (a) we have:

$M(\nu)-M(\lambda)={n\choose 2}-b(\nu^{\prime})-\left({n\choose
2}-b(\lambda^{\prime})\right)=b(\lambda^{\prime})-b(\nu^{\prime}).$

Since $\lambda_{i-1}^{*}=\lambda_{i}^{*}=2$, the domino is moved from the
second column to the first one and thus the columns sum of $\lambda$ is
decreased by $2$, while the columns sum of $\nu$ gains nothing, so that
$b(\lambda^{\prime})-b(\nu^{\prime})=2$ (see Remark 2.9).

  * •

In case (b), we have:

$b(\nu)=(i-1)(\lambda_{i-1}-1)+i(\lambda_{i}-1)+(i+1)(\lambda_{i+1}+1)+(i+2)(\lambda_{i+2}+1)+\sum_{j<i-1}{j\cdot\lambda_{j}}+\sum_{i+2<j}{j\cdot\lambda_{j}}=$

$(i-1)\lambda_{i-1}+i\lambda_{i}+(i+1)\lambda_{i+1}+(i+2)\lambda_{i+2}+4+\sum_{j<i-1}{j\cdot\lambda_{j}}+\sum_{i+2<j}{j\cdot\lambda_{j}}=b(\lambda)+4$

so that $m(\nu)-m(\lambda)=b(\nu)-b(\lambda)=4$.

∎

###  4.1. Special cases

We deal now with the case $r\leq 1$. A problem might emerge if we arrive at
the diagram $\nu=(1^{n})$ for $n$ odd, in which case, the inequality
$M(\lambda)\geq m(\nu)$ does not hold (see for example the passage from
$\lambda^{7}$ to $\lambda^{8}$ in Example 3.2). There are two possibilities of
diagrams $\lambda$ from which we might arrive at such a situation.

  1. (a)

$\lambda=(2,1^{n-2})$. This case might appear in step $(2)$ of the algorithm,
where Lemma 4.1 will not hold true since $M(\lambda)={n\choose
2}-b(\lambda^{\prime})={n\choose 2}-1$, while $m(\nu)={n\choose 2}$.
Nevertheless, this does not harm the continuity of the major index.

  2. (b)

$\lambda=(2^{2},1^{n-4})$. This can happen in the first case of step $(4)$ of
the algorithm.

In this case, we have $M(\lambda)={n\choose 2}-2$, while $m(\nu)={n\choose 2}$
(although we do have $M(\nu)-M(\lambda)=2$, the implication that
$M(\lambda)\geq m(\nu)$ which is mentioned in the parenthesized paragraph
before Fig. 9 does not hold, since the use of Lemma 2.15 with respect to $\nu$
is not allowed there).

To complete this gap, we deviate slightly from the algorithm: instead of
removing the domino in the right, we remove only the bottom square of the
domino, to get $\nu^{\prime}=(2,1^{n-2})$ in which
$M(\lambda)-M(\nu^{\prime})={n\choose 2}-\left({n\choose 2}-1\right)=1$ (note
that indeed $\nu^{\prime}\in D_{n}(1)$ since $n=2k+1$ is odd). The passage
from $\nu^{\prime}$ to $\nu=(1^{n})$ is identical to the one described in the
former case (a) (see figure 10)

$\lambda=$ &

*

$\rightarrow$ $\nu^{\prime}=$ & *

$\rightarrow$ $\nu=$

Figure 10.

We get now to the case $r=0$, here we have two types of problems.

  1. (1)

If $\lambda=(2^{2},1^{n-4})$ and $\nu=(1^{n})$, then we have
$M(\lambda)={n\choose 2}-2$, while $m(\nu)={n\choose 2}$. The value ${n\choose
2}-1$ does not appear, in accordance with the claim of Theorem 3.1. Actually,
the value ${n\choose 2}-1$ is not obtained at all due to Lemma 2.12.

  2. (2)

If $\lambda^{0}=(k,k)$ which by Lemma 2.11.1 obtains the minimum value of
major index, then since $\lambda^{0}$ is a rectangle, by Proposition 2.10,
$m(\lambda^{0})+1$ is not obtained. Again, this is in accordance with the
claim of Theorem 3.1. Note that this value will not be obtained elsewhere by
Lemma 2.12.

## Acknowledgements

The results of this paper are part of the Ph.D thesis of the second author,
supervised by Ron Adin and Yuval Roichman. We thank both of them for fruitful
conversations.

## References

  * [1] S. Billey, M. Konvalinka and J. P. Swanson, Tableau posets and the fake degrees of coinvariant algebras, Advances in Mathematics, Volume 371, 2020, 107252
  * [2] J. A. Green. The characters of the finite general linear groups, Trans. Amer. Math. Soc., 80:402–447, 1955.
  * [3] MacMahon, P. A. . The indices of permutations and the derivation therefrom of functions of a single variable associated with the permutations of any assemblage of objects, Amer. J. Math. 35 (3),(1913), 281–322.
  * [4] B. E. Sagan, The symmetric group: Representations, combinatorial algorithms, and symmetric functions, Second edition, Graduate Texts in Math., no. 203, Springer-Verlag, New York, 2001.
  * [5] R. P. Stanley, Enumerative Combinatorics, Vol. 2, Cambridge Studies in Adv. Math., no. 62. Cambridge Univ. Press, Cambridge, 1999.
  * [6] R. P. Stanley Invariants of finite groups and their applications to combinatorics, Bull. Amer. Math. Soc. (N.S.), 1(3):475–511, 1979.
  * [7] M. P. Schützenberger, La correspondence de Robinson, in: D. Foata, Ed., Combinatoire et Representation du Groupe Sym6trique, Lecture Notes in Mathematics No. 579 (Springer, Berlin, 1977) 59-113.
  * [8] R. Steinberg, A geometric approach to the representations of the full linear group over a Galois field, Trans. Amer. Math. Soc., 71:274–282, 1951.