\section{Introduction}

The inflation scenario is generally accepted for explaining the observed spatial flatness and homogeneity of the Universe.
A simple model of the scenario such as the slow-roll inflation driven by a flat inflaton potential predicts quasi de Sitter
vacuum fluctuations. These fluctuations could give rise to Gaussian and nearly scale-free metric perturbation containing both density fluctuations and gravitational waves ~\cite{olive}.

Recent cosmological observations have measured the density perturbation in an unprecedented accuracy, thus supporting the slow-roll inflation model~\cite{planckinflation}. Future observations will target at unveiling the detailed features of the metric perturbation such as non-gaussianities, the running spectral index, and the B-modes~\cite{cmb-s4}. These large-scale measurements have precisely determined the slow-roll parameters of the inflaton potential in the beginning of inflation or about 60 e-folds before the end of inflation; however, they can hardly tell us anything about inflation in a later stage or at smaller scales.

Recently, stimulated by the detection of gravitational waves from binary-black-hole mergers made by Advanced LIGO/VIRGO~\cite{LIGO}, various inflation models have been re-visited to explore the production of large curvature perturbation in a later stage of the inflationary epoch that seeds the formation of primordial black holes~\cite{pbhs,inflection}. In particular, a specific model of the single-field inflation can produce a large curvature perturbation as long as the potential has an inflection point at which the slow-roll condition is violated~\cite{inflection}. In this paper, we will scrutiny this model, concentrating on the time evolution of the superhorizon modes of the curvature perturbation.

\section{Single-field inflation model}

In a flat universe with a scale factor $a$, the equation of motion of the inflaton in a single-field inflation model with a given potential $V(\phi)$ is given by
\begin{equation}
\ddot{\phi} + 3H\dot{\phi} + \frac{dV}{d\phi} = 0,
\label{phieom}
\end{equation}
where the dot denotes the derivative with respect to the time $t$ and the Hubble parameter $H=\dot{a}/a$ is governed by the Friedmann equation,
\begin{equation}
H = \sqrt{\frac{1}{3}\left( \frac{1}{2} \dot{\phi}^2 + V \right)},
\label{heom}
\end{equation}
in which we have set the reduced Planck mass $M_p=1$. Let us introduce the Hubble flow parameters
\footnote{Originally $\epsilon$ and $\eta$ are termed as the Hubble slow-roll parameters. On cosmological scales, observations constrain them to satisfy the so-called slow-roll condition: $|\epsilon|\ll 1$ and $|\eta|\ll 1$. Here we discuss a situation with arbitrary values of $\epsilon$ and $\eta$, so for our purpose we rename them as the Hubble flow parameters.},
\begin{equation}
\epsilon = - \frac{\dot{H}}{H^2},\quad \eta = \frac{\dot\epsilon}{\epsilon H}.
\end{equation}
From Eqs.~(\ref{phieom}) and (\ref{heom}), they become
\begin{equation}
\epsilon = \frac{\dot{\phi}^2}{2 H^2},\quad \eta= -6 - \frac{2}{\dot\phi H} \frac{dV}{d\phi} + \frac{\dot\phi^2}{H^2}.
\end{equation}

In the course of inflation, quantum fluctuations of the inflaton field give rise to the curvature perturbation. In terms of the Mukhanov-Sasaki variable $v$~\cite{msvar}, the Fourier mode $v_k$ satisfies
\begin{equation}
v_k'' +\left(k^2-\frac{z''}{z}\right) v_k = 0,
\label{MS1}
\end{equation}
where the prime denotes the conformal time derivative, $d\tau=dt/a$, and $z=a\sqrt{\epsilon}$.
Here the initial condition for $v_k$ is chosen to be the standard Bunch-Davies vacuum state:
\begin{equation}
v_k\to \frac{1}{\sqrt{2k}} e^{-ik\tau},\quad\text{as}\;\frac{k}{aH} \gg 1.
\label{initialcon}
\end{equation}
We can re-write Eq.~(\ref{MS1}) in terms of $\zeta_k=v_k/(a\sqrt{2\epsilon})$ as
\begin{equation}
\zeta_k'' +\left(2 + \eta \right)aH \zeta_k' + k^2  \zeta_k = 0,
\label{MS2}
\end{equation}
and the power spectrum of the curvature perturbation is given by
\begin{equation}
\Delta^2_{\zeta_k}=\frac{k^3}{2\pi^2} \left| \zeta_k \right|^2.
\label{deltazetak}
\end{equation}

For a slow-roll inflation with $|\epsilon|\ll 1$ and $|\eta|\ll 1$,  $H$ and $\epsilon$ are approximately constant and we have $a=e^{Ht}=-1/(H\tau)$. Therefore, Eq.~(\ref{MS2}) has a simple analytic solution,
\begin{equation}
\zeta_k= \frac{\pi^{1/2}H}{2k^{3/2}}\frac{1}{\sqrt{2\epsilon}}(-k\tau)^{3/2} H^{(1)}_{3/2}(-k\tau),
\end{equation}
where only the Hankel function of the first kind appears because $\sqrt{2\epsilon}\zeta_k \to e^{-ik\tau}/(a\sqrt{2k})$ as $-k\tau\gg 1$. Further, the power spectrum of the curvature perturbation can be found by taking the late time limit ($-k\tau\ll 1$) as
\begin{equation}
\Delta^2_{\zeta_k}=\frac{H^2}{8\pi^2 \epsilon}.
\label{deltaappro}
\end{equation}

\section{A toy model with a large $\eta$}
\label{toy}

Let us consider a toy inflation model that has $|\epsilon|\ll 1$ and the temporal change of $\eta$ given by
\begin{equation}
    \eta=
    \begin{cases}
      0, & \quad\ \tau_0<\tau<\tau_1 \\
     {\bar\eta} , & \quad\ \tau_1<\tau<\tau_2 \\
      0, & \quad\ \tau_2<\tau<\tau_e
    \end{cases}
\end{equation}
where $\tau_0$ and $\tau_e$ denote the beginning and the end of inflation, respectively. Then, $H$ can be treated as a constant, $a=e^{Ht}=-1/(H\tau)$, and Eq.~(\ref{MS2}) has the solution,
\begin{equation}
   A^{-1} \zeta_k=
    \begin{cases}
      (-k\tau)^{3/2} H^{(1)}_{3/2}(-k\tau), & \quad\ \tau_0<\tau<\tau_1 \\
      (-k\tau)^\nu \left[c_1H^{(1)}_\nu(-k\tau)+c_2H^{(2)}_\nu(-k\tau)\right], & \quad\ \tau_1<\tau<\tau_2 \\
      (-k\tau)^{3/2} \left[d_1H^{(1)}_{3/2}(-k\tau)+d_2H^{(2)}_{3/2}(-k\tau)\right], & \quad\ \tau_2<\tau<\tau_e
    \end{cases}
\label{zetasoln}
\end{equation}
where
\begin{equation}
A=\pi^{1/2}H/(2k^{3/2}\sqrt{2\epsilon_0}),\quad \nu=(3+{\bar\eta})/2.
\end{equation}
The power spectrum of the curvature perturbation at the end of inflation for modes satisfying $-k\tau_e\ll 1$ is simply given by
\begin{equation}
\Delta^2_{\zeta_k}=\frac{H^2}{8\pi^2 \epsilon_0} |d_1-d_2|^2.
\label{deltad1d2}
\end{equation}

Given a $k$-mode, we can match the above piecewise functions at $\tau=\tau_1$ and $\tau=\tau_2$ to determine $c_1$, $c_2$, $d_1$, and $d_2$. First of all, it is interesting to look at a special case in which ${\bar\eta}=-6$ ($\nu=-3/2$). Using the relationships,
\begin{equation}
H^{(1)}_{-\alpha}(x)=e^{i\alpha\pi}H^{(1)}_{\alpha}(x)\quad\text{and}\quad
H^{(2)}_{-\alpha}(x)=e^{-i\alpha\pi}H^{(2)}_{\alpha}(x),
\end{equation}
we find that
\begin{equation}
c_1=i(-k\tau_1)^3,\quad c_2=0;\quad d_1=\left(\frac{\tau_2}{\tau_1}\right)^{-3},\quad d_2=0.
\label{case0}
\end{equation}
For an arbitrary ${\bar\eta}$, since we are concerned with the mode growth on superhorizon scales, we discuss three representative cases as follows.

Firstly, we consider the modes that have left the horizon before $\tau=\tau_1$ (i.e. $-k\tau_1\ll 1$).
Using the asymptotic forms at small $x$ when $\text{Re}(\alpha)>0$,
\begin{equation}
H^{(1)}_\alpha(x)\sim -H^{(2)}_{\alpha}(x)\sim - \left({i\over\pi}\right)\Gamma(\alpha) \left(x\over 2\right)^{-\alpha},
\end{equation}
we obtain that
\begin{eqnarray}
c_1-c_2\simeq\frac{\sqrt{2\pi}}{\Gamma(\nu)}2^{-\nu},
\quad d_1-d_2\simeq 1&&\quad\text{for}\quad\nu>0;\\
c_1 e^{-i\nu\pi} - c_2 e^{i\nu\pi} \simeq \frac{\sqrt{2\pi}}{\Gamma(-\nu)} 2^\nu (-k\tau_1)^{-2\nu},\quad
d_1-d_2\simeq \left(\frac{\tau_2}{\tau_1}\right)^{2\nu}&&\quad\text{for}\quad\nu<0.
\label{amplify}
\end{eqnarray}
The former case (i.e. ${\bar\eta}>-3$) means that the modes with $-k\tau_1\ll 1$ are frozen once they cross out the horizon. In the latter case, they still grow even after exiting the horizon.

Next we consider a $k$-mode such that $-k\tau_1\gg1$ and $-k\tau_2\ll 1$. Using the asymptotic forms at large $x$,
\begin{eqnarray}
H^{(1)}_{\alpha}(x)&\sim&\sqrt{\frac{2}{\pi x}}e^{i\left(x-{1\over2}\alpha\pi-{1\over4}\pi\right)},\\
H^{(2)}_{\alpha}(x)&\sim&\sqrt{\frac{2}{\pi x}}e^{-i\left(x-{1\over2}\alpha\pi-{1\over4}\pi\right)},
\end{eqnarray}
we obtain that
\begin{eqnarray}
c_1\simeq -(-k\tau_1)^{{3\over2}-\nu} e^{i\left({1\over2}\nu\pi+{1\over4}\pi\right)},\quad c_2=0;&& \label{case2}\\
d_1-d_2\simeq -\frac{2^\nu}{\sqrt{2\pi}} \Gamma(\nu) (-k\tau_1)^{{3\over2}-\nu} e^{i\left({1\over2}\nu\pi+{1\over4}\pi\right)}\quad\text{for}\quad\nu>0;&&\\
d_1-d_2\simeq  -\frac{2^{-\nu}}{\sqrt{2\pi}} \Gamma(-\nu) (-k\tau_1)^{{3\over2}+\nu} \left(\frac{\tau_2}{\tau_1}\right)^{2\nu}
e^{i\left(-{1\over2}\nu\pi+{1\over4}\pi\right)}
\quad\text{for}\quad\nu<0.&&
\end{eqnarray}

Lastly, for the $k$-mode with $-k\tau_1\gg1$ and $-k\tau_2\gg1$, we find that
\begin{eqnarray}
c_1\simeq -(-k\tau_1)^{{3\over2}-\nu} e^{i\left({1\over2}\nu\pi+{1\over4}\pi\right)},\quad c_2=0;&&\label{case3}\\
d_1\simeq \left(\frac{\tau_2}{\tau_1}\right)^{\nu-{3\over2}},\quad d_2=0.&&
\end{eqnarray}

When ${\bar\eta}=0$, in all cases we have $|d_1-d_2|^2\simeq 1$ and thus the power spectrum in Eq.~(\ref{deltad1d2}) becomes the standard slow-roll solution. When ${\bar\eta}>0$, the power spectrum is suppressed. When ${\bar\eta}<0$, due to the mode growth on superhorizon scales, the power spectrum is amplified, for example in Eq.~(\ref{amplify}), by a factor,
\begin{equation}
\left(\frac{\tau_2}{\tau_1}\right)^{4\nu} = e^{-4\nu\Delta N},
\label{amplifyN}
\end{equation}
where the e-folding is defined by $N=\int_0^t H(t')dt'$ and $\Delta N=N_2-N_1$ is the change of e-folds from $\tau_1$ to $\tau_2$.

\section{Two-point correlation function}

One may ask whether the superhorizon growth violates the causality. To answer this question, we examine the correlation length by computing the equal-time two-point correlation function. In order to save the causality, the correlation length should always be smaller than the causal horizon.

The two-point correlation function is defined by
\begin{equation}
\langle\zeta({\vec x},\tau) \zeta({\vec x}',\tau') \rangle
= \int \frac{d^3{\vec k}}{(2\pi)^3} \zeta_k(\tau) \zeta_k^*(\tau') e^{i{\vec k}\cdot({\vec x}-{\vec x}')},
\label{2ptfct}
\end{equation}
where $\zeta_k$ is the solution in Eq.~(\ref{zetasoln}).  It suffices to consider the function in the period $\tau_1<\tau<\tau_2$ when ${\bar\eta}\neq 0$. The results in the previous section suggest that we can take $c_1$ of the form $c_1= Ck^\beta$ and set $c_2=0$. Then, Eq.~(\ref{2ptfct}) becomes
\begin{eqnarray}
&&\langle\zeta({\vec x},\tau) \zeta({\vec x}',\tau') \rangle \nonumber \\
&=&  \frac{\pi H^2}{8\epsilon_0}  |C|^2 (\tau\tau')^{\nu}
\int \frac{d^3{\vec k}}{(2\pi)^3} k^{2(\beta+\nu)-3} H^{(1)}_\nu(-k\tau) H^{(1)*}_\nu(-k\tau')
e^{i{\vec k}\cdot({\vec x}-{\vec x}')} \nonumber \\
&=&  \frac{H^2}{4\pi^3\epsilon_0}  |C|^2 \frac{(\tau\tau')^{\nu}}{\left|{\vec x}-{\vec x}'\right|}
\int_0^\infty dk\, k^{2(\beta+\nu)-2} \sin\left(k\left|{\vec x}-{\vec x}'\right|\right)
K_\nu(ik\tau) K_\nu(-ik\tau'). \nonumber
\end{eqnarray}
In most of the cases, as found in Eqs.~(\ref{case0}), (\ref{case2}), and (\ref{case3}),
we have $|C|^2=|\tau_1|^{3-2\nu}$ and $\beta=3/2-\nu$. Hence, we obtain that
\begin{eqnarray}
&&\langle\zeta({\vec x},\tau) \zeta({\vec x}',\tau') \rangle \nonumber \\
&=&  \frac{H^2}{32\pi^2\epsilon_0} \left(\frac{\tau\tau'}{\tau_1^2}\right)^{\nu-{3\over2}}
\Gamma\left({3\over2}+\nu\right) \Gamma\left({3\over2}-\nu\right)
{_2}F_1\left({3\over2}+\nu,{3\over2}-\nu;2;1-{1\over w}\right),\label{hyper}
\end{eqnarray}
where $|\nu|<3/2$ and
\begin{equation}
w=\frac{4\tau\tau'}{-(\tau-\tau')^2+\left|{\vec x}-{\vec x}'\right|^2}\,.
\nonumber
\end{equation}
For any two spatial points separated by a distance much greater than the horizon, i.e. $\left|{\vec x}-{\vec x}'\right|\gg |\tau|$, we can set $\tau=\tau'$ and take the limit $w\rightarrow 0^+$ in Eq.~(\ref{hyper}) to obtain the equal-time correlation,
\begin{equation}
\langle\zeta({\vec x},\tau) \zeta({\vec x}',\tau) \rangle =
\frac{H^2}{8\pi^{5/2}\epsilon_0} \left(\frac{\tau^2}{\tau_1^2}\right)^{\nu-{3\over2}}
\left(\frac{\tau^2}{\left|{\vec x}-{\vec x}'\right|^2}\right)^{{3\over2}-|\nu|}
 \Gamma\left(|\nu|\right) \Gamma\left({3\over2}-|\nu|\right),
\end{equation}
which diminishes with the two-point separation when $|\nu|<3/2$ ($0>{\bar\eta}>-6$).
This means that although the modes are growing on superhorizon scales,
the growth is causal in the sense that the correlation length is safely smaller than the horizon.
It is worth noting that by means of the toy model, we have shown that the causality condition requires that
${\bar\eta} \gtrsim -6$. In the next section, we will use a realistic inflation model to further examine this lower limit.

\begin{figure}[htp]
\centering
\includegraphics[width=0.8\textwidth]{fig1}
\caption{The inflaton potential $V(\phi)$ in Eq.~(\ref{infpot}), normalized by $V_0$. All dynamical variables in this figure and in the following figures are rescaled by the reduced Planck mass, $M_p=2.435\times 10^{18}$ GeV. }
\label{fig1}
\end{figure}

\begin{figure}[htp]
\centering
\includegraphics[width=0.8\textwidth]{fig2}
\caption{Evolution of $H$ (dashed line) and $\phi$ (solid line) against e-folds $N$, with
$H_0=7.65\times10^{-8}$, $\phi_0=8.48\times 10^{-4}$, and ${\dot\phi}_0=1.90\times 10^{-11}$.
Note that the inflation ends at $N\sim 60$.}
\label{fig2}
\end{figure}

\begin{figure}[htp]
\centering
\includegraphics[width=0.8\textwidth]{fig3}
\caption{Evolution of $\epsilon = - \dot{H}/H^2$ (dashed line) and $\eta = \dot\epsilon/(\epsilon H)$ (solid line).}
\label{fig3}
\end{figure}

\begin{figure}[htp]
\centering
\includegraphics[width=0.8\textwidth]{fig4}
\caption{Power spectrum of the curvature perturbation. The thick solid curve is the exact numerical result from solving Eq.~(\ref{MS2}). The dashed curve is an approximate result obtained by Eq.~(\ref{deltaappro}). The thin solid curve is the numerical result with $\zeta_k$ evaluated at the horizon-crossing time for each $k$-mode.}
\label{fig4}
\end{figure}

\begin{figure}[htp]
\centering
\includegraphics[width=0.8\textwidth]{fig5}
\caption{Snapshots of the evolution of the curvature perturbation in the course of the inflation, timed by successive e-folds $N=32,33,34,35,36,37,38,60$. When drawing a curve for the snapshot at time $N$, the $k$-modes that have not yet crossed out the horizon are omitted, so the curve ends at the e-folding $N$. }
\label{fig5}
\end{figure}

\section{Inflation with an inflection point}

In order to have a large value of $\eta$, we utilize an inflation model with an inflection point. The potential is given by
\begin{equation}
V(\phi)=V_0 \left[ 1+ c_1 \frac{\phi}{c_\Lambda} + \frac{c_3}{3!} \left(\frac{\phi}{c_\Lambda}\right)^3  + \frac{c_4}{4!} \left(\frac{\phi}{c_\Lambda}\right)^4 + \frac{c_5}{5!} \left(\frac{\phi}{c_\Lambda}\right)^5 \right].
\label{infpot}
\end{equation}
We tune the above potential parameters and the initial conditions of the inflaton motion, given by
\begin{eqnarray}
&&c_\Lambda = 0.3,\quad V_0 = 1.7574 \times 10^{-14}, \nonumber \\
&&c_1 = -0.7278 \times 10^{-4},\; c_3 = -0.52,\; c_4 = 1.0,\; c_5 = -0.6407; \nonumber \\
&&\phi_0 = 8.5 \times 10^{-4},\quad {\dot\phi}_0 = 1.9 \times 10^{-11};
\end{eqnarray}
such that at near the inflection point the inflaton undergoes an ultra-slow-roll motion, while the inflaton acceleration is non-negligible. In addition, the amplitude of the power spectrum of the curvature perturbation as well as the levels of the scalar spectral index, the tensor-to-scalar ratio, and the scalar spectral index running on large cosmological scales are all consistent with the Planck measurements:
$\Delta_{\zeta_k}^2\simeq 2.4\times 10^{-9}$, $n_s\simeq 0.97$, $r<0.06$, and $|dn_s/d\ln k|<0.013$~\cite{planckinflation,bk18}.

Figure~\ref{fig1} is the inflaton potential with an inflection point. Figure~\ref{fig2} depicts the evolution of the Hubble parameter and the inflaton mean field. We can see that the inflaton field is almost stationary near the inflection point at $\phi\simeq 0.95$ and hence prolongs the duration of inflation by about $25$ e-folds. The evolution of the Hubble flow parameters is shown in Fig.~\ref{fig3}, where a sudden drop in $\epsilon$ at $N\sim 32$ results in a large decrease of $\eta$. The value of $\eta$ has reached about $-7$ at $N\sim 35$, which is close to the predicted lowest limit of $\eta$ by the causality argument in the previous section.

We numerically solve for $\zeta_k$ in Eq.~(\ref{MS2}) under the initial condition~(\ref{initialcon}). The power spectrum of the curvature perturbation is then given by Eq.~(\ref{deltazetak}) with $\zeta_k$ evaluated at the end of inflation. In the present consideration, we specifically compare the exact numerical result with that directly obtained by the approximate formula in Eq.~(\ref{deltaappro}). Figure~\ref{fig4} shows that the approximate formula underestimates the rapid growth of the curvature perturbation when $\eta$ has a large negative value. Also, in order to highlight the superhorizon growth, we have produced a power spectrum using Eq.~(\ref{deltazetak}) in which $\zeta_k$ takes the value when the $k$-mode crosses the horizon. This power spectrum, denoted by the thin solid curve in Fig.~\ref{fig4}, matches the one from the approximate formula. In Fig.~\ref{fig5}, we show the growth of the curvature perturbation in the course of the inflation.

In Sec.~\ref{toy}, based on the toy model, we have shown that the growth of the superhorizon modes can boost up the power spectrum by a factor of $e^{-4\nu\Delta N}$. We can check this with the numerical results as follows. In Fig.~\ref{fig3}, the evolution profile of $\eta$ can be roughly captured by the shape of a square well with a depth of about $-7$ at $N\sim 35$ and a width of $\Delta N \simeq 3$. This gives rise to a boost factor of about $e^{24}\sim10^{12}$, which matches fairly well to the large power increase from the thin solid line to the exact power spectrum at $N\sim 34$ in Fig.~\ref{fig4}. Figure~\ref{fig5} depicts how the growth of the superhorizon modes during the short interval of a few e-folds at $N\sim 34$ renders the large boost to the power spectrum.

\section{Conclusions}

We have analyzed the evolution of the superhorizon modes of the curvature perturbation during single-field inflation in a period when the slow-roll condition is violated. The conservation of the superhorizon curvature perturbation, derived in the assumption of the slow-roll condition~\cite{conservation}, is no longer valid, as the superhorizon modes can grow rapidly in an ultra-slow-roll phase during which the speed of the inflaton is extremely low, while the acceleration is large. Nevertheless, the principle of causality implies an upper bound on the growth rate.

We have illustrated the superhorizon growth in a full numerical calculation with the inflaton potential having an inflection point. The inflation undergoes an ultra-slow-roll phase at the inflection point such that the value of $\epsilon$ is extremely small whereas $\eta$ can reach a negative value as low as about $-7$. This value is consistent with the lower limit of $\eta\sim -6$ estimated by an analytic study in the present work. As a consequence, the curvature perturbation can grow enormously to a peak value of $\Delta^2_{\zeta_k}\sim 0.01$ at $N\sim 34$, which serves as seed for forming primordial black holes of mass about $10^{-8} M_\odot$.

Our analytic study can be readily used to assess the growth of the curvature perturbation during ultra-slow-roll inflation for estimating the abundances and the mass ranges of primordial black holes formed when the curvature perturbation re-enters the horizon.


\section{Introduction}
\label{sec.intro}

Among lattice models of strongly correlated fermions, the $t$-$J$ model holds a prominent position
due to its role in the theoretical description of high-temperature superconductivity and as an
example of a simple model featuring spin-charge separation~\cite{Sc87,ZR88,EKS92}. In this model
each lattice site can be either empty or occupied by one fermion, which interacts with its nearest
neighbors through spin exchange and charge repulsion and can also hop between contiguous lattice
sites.. In the one-dimensional case the $t$-$J$ model is of particular interest, as it is both
supersymmetric and exactly solvable through the nested Bethe ansatz when its parameters are
suitably related~\cite{Sc87,Su75,BB90,KY90,Sa91,EK92}.
In the early nineties, Kuramoto and Yokoyama~\cite{KY91,Ku93} introduced a long-range version of
the supersymmetric $t$-$J$ model featuring $1/r^2$ interactions, which reduces to the $\mathrm{su}(2)$
Haldane--Shastry (HS) chain \cite{Ha88,Sh88} in the high density limit (i.e., when all the sites
are occupied). Among other interesting features, the Kuramoto--Yokoyama (KY) model exhibits strong
spin-charge separation, in the sense that at low temperatures the spin and charge velocities are
respectively independent of the charge density and the magnetization.
At low energies the KY model is known to be a Luttinger liquid~\cite{KY91b}, with spin and charge
excitations independently described by a $c=1$ conformal field theory (CFT).

The supersymmetric KY model is exactly solvable through the asymptotic Bethe ansatz pioneered by
Sutherland and Shastry~\cite{SS93}, as its energies can be found in principle by solving a system
of transcendental equations in the asymptotic momenta~\cite{Ka92,Ka92b}. However, this method does
not completely determine the spectrum, since it does not provide complete information on the
degeneracy of each level. Based on numerical calculations, Wang et al.~\cite{WLC92} proposed an
empirical description of the degeneracies of the spectrum of the $\mathrm{su}(2)$ KY model reminiscent of
the rule for filling the border strips associated to Haldane's motifs~\cite{HHTBP92,Ha93,BBHS07}.
This description, however, is known to be incorrect in certain situations, although the needed
corrections vanish in the thermodynamic limit~\cite{WLC92}. Inspired by the equivalence between
the $\mathrm{su}(2)$ KY model and the $\mathrm{su}(1|2)$ supersymmetric HS chain (up to a term proportional to the
total electric charge), Saiga and Kuramoto~\cite{SK99} conjectured a description of the former
model's spectrum essentially in terms of $\mathrm{su}(1|2)$-supersymmetric Haldane motifs, which accounted
for the numerical results for $N\le 16$ spins. To the best of our knowledge, this conjecture has
remained unproved in the literature.

In this paper we address and completely solve the problem of finding a full description of the
spectrum of the supersymmetric $\mathrm{su}(m)$ KY model with a general chemical potential term for a
finite number of sites, including the determination of the levels' degeneracies and spin content.
In particular, our results provide a rigorous proof of the Saiga--Kuramoto conjecture for
arbitrary $m$ and $N$. Our approach, which we shall now briefly summarize, is new and bypasses the
usual machinery of the asymptotic Bethe ansatz, transfer matrix, Yangian highest-weight states,
etc. Indeed, we start by establishing the precise connection between the $\mathrm{su}(m)$ KY model and the
$\mathrm{su}(1|m)$ supersymmetric HS spin chain with a suitable chemical potential term for arbitrary $m$,
thus generalizing the well-known result for $m=2$. Since the partition function of the latter
chain was recently evaluated in Ref.~\cite{FGLR18} by using Polychronakos's freezing trick, this
immediately yields a novel closed formula for the partition function of the $\mathrm{su}(m)$ KY model. A
remarkable property of this partition function is that it can be recast as the partition function
of an equivalent (inhomogeneous) vertex model with simple interactions~\cite{BBH10}. We show that
both partition functions are polynomials in appropriate variables, whose coefficients are nothing
but the corresponding restricted partition functions on subspaces with well-defined magnon
content. This crucial observation provides a closed formula for the restricted partition function
of the $\mathrm{su}(m)$ KY model (in the presence of an external magnetic field and a chemical potential
term) on each of these subspaces. Finally, by analyzing the restricted partition function of the
equivalent vertex model we are able to express the spectrum of the $\mathrm{su}(m)$ KY model in each
subspace with well-defined magnon content in terms of suitably restricted supersymmetric Haldane
motifs and their corresponding Young tableaux.
This yields a complete and rigorous description of the spectrum in the latter subspaces, including
a systematic way for determining the degeneracy of each level, which implies the Saiga--Kuramoto
conjecture as a particular case.

It should be noted that, while the traditional freezing trick allows one to compute the partition
function (and, in principle, the spectrum) of certain integrable systems, it does not provide any
information about the corresponding wave functions. On the other hand, the analysis of the
spectrum of the $\mathrm{su}(m)$ KY model by the method described above, extends the freezing trick and,
more importantly, the equivalence to a vertex model, to subspaces of the Hilbert space with
well-defined magnon content. It is in fact this connection with a restricted vertex model what
makes it possible to identify all energy eigenvalues within any subspace of the Hilbert space with
well-defined magnon content. In other words, our approach not only yields the complete spectrum of
the $\mathrm{su}(m)$ KY model but also the magnon numbers or spin content of the corresponding wave
functions.

Our approach yields several additional results that we shall now briefly discuss. In the first
place, as a consequence of the general discussion of the equivalence of the $\mathrm{su}(m)$ KY model to
an $\mathrm{su}(1|m)$ supersymmetric HS chain with a suitable chemical potential term, we construct a new
model of KY type with general elliptic interactions which can be mapped to a corresponding
$\mathrm{su}(1|m)$ generalization of Inozemtsev's chain~\cite{In90}. This model certainly deserves further
study, since it smoothly interpolates between the standard (nearest-neighbors) $t$-$J$ model and
the (long-range) KY model. Secondly, as an application of the description of the spectrum of the
KY model in terms of supersymmetric Young tableaux, we determine the ground state of the spin
$1/2$ model in an external magnetic field in the thermodynamic limit. In particular, we give a
complete description of the different ground state phases, characterized by their spin content
---i.e., $\mathrm{su}(1|2)$, $\mathrm{su}(1|1)$ and $\mathrm{su}(0|2)$, apart from the trivial phases consisting of only
holes or fermions of one type--- in terms of the magnetic field strength and the charge chemical
potential. This description goes beyond previously known results, which are restricted to the
genuinely~$\mathrm{su}(1|2)$ phase. In particular, we show that the strong spin-charge separation
characteristic of the long-range $t$-$J$ model at low temperatures~\cite{KK95} occurs in all
nontrivial phases.

The paper is organized as follows. In Section~\ref{sec.models} we introduce the model and show its
equivalence to a supersymmetric HS chain with suitable chemical potential terms. We also introduce
the model's elliptic version and discuss its connection with the supersymmetric extension of
Inozmetsev's elliptic chain. The model's partition function is then computed in
Section~\ref{sec.PF} by exploiting its equivalence to a supersymmetric HS chain. In
Section~\ref{sec.motifs} we recast the partition function as that of a suitable inhomogeneous
vertex model and derive the model's restricted partition function on subspaces with well-defined
magnon content. As explained above, this yields an explicit and complete description of the
spectrum on the latter subspaces, including each level's degeneracy, in terms of suitable
supersymmetric Young tableaux. Section~\ref{sec.GS} is devoted to a complete analysis of the
ground state phases of the spin~$1/2$ model in an external magnetic field in the thermodynamic
limit. In Section~\ref{sec.conc} we present our conclusions and outline several future
developments. The paper ends with a short technical appendix in which we present the proof of a
new result regarding the degeneracy of reverse motifs.

\section{The models}\label{sec.models}

We shall deal in this paper with a class of~$\mathrm{su}(m)$ $t$-$J$ type models, consisting of a
one-dimensional lattice with $N$ sites each of which can be either empty or occupied by a single
fermion with $m$ internal degrees of freedom. We shall be mainly interested in long-ranged models
of the latter type, in which the spin and charge interactions among the fermions and their hopping
amplitude involve all possible pairs of lattice sites. More precisely, we shall take as the
model's Hamiltonian
\begin{equation}
  \fl
  H_0=\sum_{\mathpalette\mathclapinternal{1\le i<j\le N}}{\mathcal P}\bigg\{-t_{ij}\sum_{\sigma=1}^m(c^\dagger_{i\sigma}c_{j\sigma}
  +c^\dagger_{j\sigma}c_{i\sigma})+2J_{ij}\mathbf{T}_i\cdot\mathbf{T}_j-2V_{ij}n_in_j\bigg\}{\mathcal P}\equiv
  \sum_{\mathpalette\mathclapinternal{1\le i<j\le N}}H_{ij}\,,
  \label{tJgen}
\end{equation}
or equivalently\footnote{Here and in what follows, unless otherwise stated sums and products over
  Latin indexes run over the set~$1,\dots,N$, while Greek indices range from~$1$ to~$m$.}
\begin{equation}
  \label{tJgen2}
  H_0=\sum_{i\ne j}{\mathcal P}\bigg\{-t_{ij}\sum_\sigmac^\dagger_{i\sigma}c_{j\sigma}+
  J_{ij}\mathbf{T}_i\cdot\mathbf{T}_j-V_{ij}n_in_j\bigg\}{\mathcal P}\,,
\end{equation}
where~$t_{ij}=t_{ji}$, $J_{ij}=J_{ji}$, $V_{ij}=V_{ji}$ are real constants. We shall also assume
that the model~\eqref{tJgen}-\eqref{tJgen2} is \emph{translation-invariant}, i.e., that
\begin{equation}\label{tJV}
  t_{ij}=t(i-j)\,,\quad
  J_{ij}=J(i-j)\,,\quad V_{ij}=V(i-j)
\end{equation}
with
\begin{equation}\label{trinv}
  t(x)=t(-x)=t(N-x)\,,
\end{equation}
and similarly for $J(x)$, $V(x)$. In the latter equations~$c^\dagger_{i\sigma}$ (respectively~$c_{i\sigma}$)
denotes the operator creating (resp.~destroying) a fermion of type~$\sigma$ at site~$i$ and
$n_i=\sum_\sigma n_{i\sigma}$, where $n_{i\sigma}=c^\dagger_{i\sigma}c_{i\sigma}$, is the total number of fermions at
site~$i$. The operator ${\mathcal P}$ is the projector onto single-occupancy states, in which each site is
occupied by at most one fermion. Finally, $\mathbf{T}_i\equiv(T_i^1,\dots,T_i^{m^2-1})$, where~$T_i^r$ is
the~$r$-th $\mathrm{su}(m)$ Hermitian generator in the fundamental representation acting on the $i$-th
site (with a suitable normalization that we shall specify below). Thus the first term
(proportional to $t_{ij}$) in Eqs.~\eqref{tJgen}-\eqref{tJgen2} accounts for the hopping of
fermions between sites~$i$ and~$j$, while the last two terms respectively model the spin
(exchange) and charge interaction between the latter sites.

The Hamiltonian~\eqref{tJgen}-\eqref{tJgen2} encompasses several well-known models, which we shall
briefly review. To begin with, note that a nearest-neighbors version of the
Hamiltonian~\eqref{tJgen}-\eqref{tJgen2} is obtained by taking~$t(x)$, $J(x)$, $V(x)$ proportional
to the $N$-periodic extension of
\begin{equation}\label{tde}
  \delta_{1,x}+\delta_{N-1,x}\,,\qquad 1\le x\le N-1\,.
\end{equation}
When~$m=2$, this is the original $t$-$J$ model introduced in Ref.~\cite{Sc87}. On the other hand,
the long-ranged supersymmetric Kuramoto--Yokoyama model~\cite{KY91,Ku93} follows from
Eqs.~\eqref{tJgen2}-\eqref{tJV} when
\[
  t(x)=J(x)=4V(x)=\frac{t\pi^2}{N^2}\sin^{-2}\bigl(\pi x/N\bigr)\,,
\]
where~$t$ is a positive real parameter. More generally, when~$m=2$ the
model~\eqref{tJgen2}-\eqref{tJV} is supersymmetric provided that $t(x)=J(x)=4V(x)$. In fact, one
of our aims in this section is to generalize the latter result to the $\mathrm{su}(m)$ case with arbitrary
$m>2$. As we shall next see, the key idea in this respect is to map the original
Hamiltonian~\eqref{tJgen} to that of a suitable supersymmetric spin chain in which the holes are
regarded as bosons.

More precisely, consider a one-dimensional lattice (\emph{spin chain}) each of whose sites are
occupied either by a boson or an~$\mathrm{su}(m)$ fermion. Thus the model's Hilbert space
is~$\hat{\mathcal H}=\otimes_{i=1}^N\hat{\mathcal H}_i$, where~$\hat{\mathcal H}_i$ is the linear span of the one-particle
states $b^\dagger_i\ket{\hat\Omega}_i$, $f^\dagger_{i\sigma}\ket{\hat\Omega}_i$ with $\sigma=1,\dots,m$, $b^\dagger_i$,
$f^\dagger_{i\sigma}$ are the operators creating respectively a boson and a fermion of type~$\sigma$ at the
$i$-th site, and~$\ket{\hat\Omega}_i$ is the vacuum. Similarly, denote by ${\mathcal H}=\otimes_{i=1}^N{\mathcal H}_i$
the Hilbert space of the original model~\eqref{tJgen}, where~${\mathcal H}_i$ is spanned by the
vacuum~$\ket\Omega_i$ and the one-particle states~$c^\dagger_{i\sigma}\ket{\Omega}_i$. We now introduce the
unitary mapping~$\varphi:{\mathcal H}\to\hat{\mathcal H}$ defined by
\[
  \varphi\ket\Omega_i=b^\dagger_i\ket{\hat\Omega}_i\,,\quad
  \varphi\big(c^\dagger_{i\sigma}\ket\Omega_i\big)=f^\dagger_{i\sigma}\ket{\hat\Omega}_i\,.
\]
This mapping induces a natural way of associating to each linear operator $A$ acting on~${\mathcal H}$ a
corresponding linear operator~$\hat A=\varphi A\varphi^{-1}=\varphi A \varphi^\dagger$ acting on $\hat{\mathcal H}$. Note,
in particular, that
$(A^\dagger)^{\textstyle\hat{}}=\varphi A^\dagger \varphi^\dagger\equiv\hat A^\dagger$. It is
straightforward to check that under this correspondence $\hat c_{i\sigma}=b^\dagger_if_{i\sigma}$, since both
operators agree on the canonical basis of~$\hat{\mathcal H}_i$. Note than on~${\mathcal H}_i$ we
have~$c_{i\sigma}={\mathcal P} c_{i\sigma}{\mathcal P}$, so that we can also write
\begin{equation}\label{PcP}
  ({\mathcal P} c_{i\sigma}{\mathcal P})^{\textstyle\hat{}}=b^\dagger_if_{i\sigma}\equiv X_i^{0\sigma},
\end{equation}
and therefore (since ${\mathcal P}$ is Hermitian, being a projector)
\begin{equation}\label{PdcP}
  ({\mathcal P}c^\dagger_{i\sigma}{\mathcal P})^{\textstyle\hat{}}=\hat c_{i\sigma}^\dagger=f^\dagger_{i\sigma}b_i\equiv X_i^{\si0}\,.
\end{equation}
We shall also need in the sequel the relations
  \begin{eqnarray}
    \smash{\big({\mathcal P} c^\dagger_{i\sigma}c_{i\sigma'}{\mathcal P}\big)}^{\textstyle\hat{}}
    &=f^\dagger_{i\sigma}f_{i\sigma'}\equiv
      X_i^{\sigma\si'}\,,
      \label{Xissp}\\
    {[{\mathcal P}(1-n_i){\mathcal P}]}^{\,\textstyle\hat{}}&=b^\dagger_ib_i\equiv X_i^{00}\,,
                                          \label{niXi00}
  \end{eqnarray}
which easily follow from the previous ones. For instance, taking into account
that~$c_{i\sigma'}{\mathcal P}={\mathcal P} c_{i\sigma'}{\mathcal P}$ and~${\mathcal P}^2={\mathcal P}$ we have
\begin{eqnarray*}
  \smash{\big({\mathcal P} c^\dagger_{i\sigma}c_{i\sigma'}{\mathcal P}\big)}^{\textstyle\hat{}}
  &=\smash{\big({\mathcal P} c^\dagger_{i\sigma}{\mathcal P}\cdot{\mathcal P} c_{i\sigma'}{\mathcal P}\big)}^{\textstyle\hat{}}
  \\
  &=f^\dagger_{i\sigma}b_ib^\dagger_if_{i\sigma'}=X_i^{\sigma\si'}+f^\dagger_{i\sigma}b^\dagger_if_{i\sigma}b_i
    =X_i^{\sigma\si'},
\end{eqnarray*}
since~$f_{i\sigma}b_i=0$ on~$\hat{\mathcal H}_i$.

Consider next the $\mathrm{su}(1|m)$ supersymmetric permutation
operators~$P_{ij}^{(1|m)}:\hat{\mathcal H}\to\hat{\mathcal H}$ (with~$i<j$), whose action on the canonical basis
\[
  \ket{\sigma_1\cdots \sigma_N}\equiv a^\dagger_{1\sigma_1}\cdotsa^\dagger_{N\sigma_N}\ket{\hat\Omega}\,,\quad
  \sigma_1,\dots,\sigma_N\in\{0,\dots,m\},
\]
where~$a_{i0}=b_i$ and~$a_{i\sigma}=f_{i\sigma}$ for~$\sigma\ge1$
and~$\ket{\hat\Omega}=\otimes_i\ket{\hat\Omega}_i$ is the global vacuum, is given by~\cite{BB06}
\[
  P_{ij}^{(1|m)}\ket{\cdots\sigma_i\cdots\sigma_j\cdots}=\epsilon({\boldsymbol{\si}})\ket{\cdots\sigma_j\cdots\sigma_i\cdots}\,.
\]
The sign~$\epsilon({\boldsymbol{\si}})$ is $1$ (respectively $-1$) if $\sigma_i=\sigma_j=0$ (resp.~$\sigma_i,\sigma_j\ge1$),
while for $\sigma_i\sigma_j=0$ and $\sigma_i\ne\sigma_j$ it is equal to the number of fermionic spins $\sigma_k$
with~$i+1\le k\le j-1$. It is well known~\cite{Ha93,HB00,BB06} that~$P_{ij}^{(1|m)}$ can be
expressed in terms of creation and annihilation operators as follows:
\[
  P_{ij}^{(1|m)}=\sum_{\alpha,\beta=0}^ma^\dagger_{i\alpha}a^\dagger_{j\beta}a_{i\beta}a_{j\alpha}=
  \sum_{\alpha,\beta=0}^m(-1)^{p(\beta)}X_{i}^{\alpha\beta}X_j^{\beta\alpha}\,,
\]
where~$p(0)=0$ and $p(\sigma)=1$ for~$\sigma\ge1$. We thus have
\begin{equation}\label{Pij1m}
  P_{ij}^{(1|m)}=X_i^{00}X_j^{00}+\sum_\sigma(X_i^{\si0}X_j^{0\sigma}-X_i^{0\sigma}X_j^{\si0})
  -P_{ij}\,,
\end{equation}
where
\[
  P_{ij}=\sum_{\sigma,\sigma'}X_{i}^{\sigma\si'}X_j^{\sigma'\!\sigma}
\]
is the ordinary permutation operator when acting on purely fermionic states.

Our next goal is to relate the product~$\mathbf{T}_i\cdot\mathbf{T}_j$ appearing in the
Hamiltonian~\eqref{tJgen} with the supersymmetric permutation operator~$P_{ij}^{(1|m)}$. To this
end, note first of all that the components~$T_i^r$ of~$\mathbf{T}_i$ are defined in the usual way as
\begin{equation}\label{Tir}
  T_i^r=\sum_{\sigma,\sigma'}T^r_{\sigma\si'}c^\dagger_{i\sigma}c_{i\sigma'}\,,
\end{equation}
where the complex numbers~$T^r_{\sigma\si'}$ are the matrix elements of the~$r$-th (Hermitian)
generator of~$\mathrm{su}(m)$ in the fundamental representation. We shall normalize the~$m\times m$
matrices~$T^r\equiv \big(T^r_{\sigma\si'}\big)_{1\le\sigma,\sigma'\le m}$ so that
\[
  \operatorname{tr}(T^rT^s)=\frac12\,\delta_{rs}\,.
\]
In particular, when~$m=2$ the operator $T^r$ can be taken as the usual spin~$1/2$
operator~$S^r=\sigma^r/2$, where~$\sigma^r$ is the $r$-th Pauli matrix. In order to
relate~$\mathbf{T}_i\cdot\mathbf{T}_j$ with~$P_{ij}^{(1|m)}$ we shall make use of the identity
\begin{equation}\label{comp}
  2\sum_{r=1}^{m^2-1}(T^r)_{\sigma\si'}(T^r)_{\mu\mu'}=\delta_{\sigma\mu'}\delta_{\sigma'\mu}
  -\frac1m\,\delta_{\sigma\si'}\delta_{\mu\mu'}\,,
\end{equation}
which is a direct consequence of the completeness of the generators $T^r$ together with the
identity matrix. From Eqs.~\eqref{Xissp}, \eqref{Tir} and~\eqref{comp} we obtain
\begin{eqnarray}
  \fl
  2({\mathcal P}\mathbf{T}_i\cdot\mathbf{T}_j{\mathcal P})^{\textstyle\hat{}}
  &=2\sum_{r=1}^{m^2-1}\sum_{\substack{\sigma,\sigma'\cr\mu,\mu'}}T^r_{\sigma\si'}T^r_{\mu\mu'}
    X_i^{\sigma\si'}
    X_j^{\mu\mu'}=\sum_{\sigma,\sigma'}X_i^{\sigma\si'}X_j^{\sigma'\!\sigma}-\frac1m\sum_{\sigma,\mu}X_i^{\sigma\si}X_i^{\mu\mu}
    \nonumber\\
  \fl
  &= P_{ij}-\frac1m\,\hat n_i\hat n_j\,,
    \label{TiTjPij}
\end{eqnarray}
where~$\hat n_k=\sum_\sigma f^\dagger_{k\sigma}f_{k\sigma}$ denotes the total number of fermions (created by the
operators $f^\dagger_{k\sigma}$) at the~$k$-th site. From Eq.~\eqref{Pij1m} for the supersymmetric
permutation operator~$P_{ij}^{(1|m)}$ we obtain, after some algebra,
\begin{equation}\label{PijAij}
  P_{ij}^{(1|m)}+\hat n_i+\hat n_j-1=\hat A_{ij}\,,
\end{equation}
where
\[
  A_{ij}={\mathcal P}\bigg[\sum_\sigma(c^\dagger_{i\sigma}c_{j\sigma}+c^\dagger_{j\sigma}c_{i\sigma})-2\mathbf{T}_i\cdot\mathbf{T}_j
  +\bigg(1-\frac1m\bigg)n_in_j\bigg]{\mathcal P}.
\]
Comparing with Eq.~\eqref{tJgen} we deduce that~$H_{ij}$ will be proportional to~$A_{ij}$ provided
that
\begin{equation}
  \label{SUSY}
  t(x)=J(x)=2\bigg(1-\frac1m\bigg)^{\!\!-1}V(x)\,,
\end{equation}
and in that case
\[
  \hat H_{ij}=t_{ij}\big(1-P_{ij}^{(1|m)}-\hat n_i-\hat n_j\big)\,.
\]
In other words, when condition~\eqref{SUSY} is satisfied, i.e., when~$H_0$ is of the form
\begin{equation}
  H_0=\sum_{i<j}t_{ij}{\mathcal P}\bigg[-\sum_{\sigma}(c^\dagger_{i\sigma}c_{j\sigma}+c^\dagger_{j\sigma}c_{i\sigma})
  +2\mathbf{T}_i\cdot\mathbf{T}_j-\left(1-\tfrac1m\right)n_in_j\bigg]{\mathcal P}\,,
  \label{tJSUSY}
\end{equation}
the corresponding Hamiltonian~$\hat H_0$ is given by
\[
  \hat H_0=\sum_{i<j}t_{ij}(1-P_{ij}^{(1|m)}-\hat n_i-\hat n_j)\\
  =\sum_{i<j}t_{ij}(1-P_{ij}^{(1|m)})-\sum_{i\ne j}t_{ij}\hat n_j\,.
\]
Note that so far we have not used the translation-invariance conditions~\eqref{tJV}-\eqref{trinv},
so that the previous result is valid in full generality. On the other hand, when the
model~\eqref{tJgen} is translation-invariant we can use Eqs.~\eqref{tJV}-\eqref{trinv} to further
simplify the last term in the previous equation. Indeed, in this case
\[
  \sum_{i\ne j}t_{ij}\hat n_j=\sum_j\bigg(\sum_{i;i\ne j}t_{ij}\bigg)\hat n_j\,,
\]
with
\[
  \fl
  \sum_{i;i\ne j}t_{ij}=\sum_{k=1-j}^{-1}t(k)+\sum_{k=1}^{N-j}t(k)
  = \sum_{k=N-j+1}^{N-1}t(k-N)+\sum_{k=1}^{N-j}t(k) = \sum_{k=1}^{N-1}t(k)\equiv t_0\,,
\]
so that
\begin{equation}
  \hat H_
  =\sum_{i<j}t_{ij}(1-P_{ij}^{(1|m)})-t_0{\mathcal F},
  \label{hatH}
\end{equation}
where
\[
  {\mathcal F}\equiv\sum_{i}\hat n_{i}
\]
is the total number of fermions.

Summarizing, we have shown that the translation-invariant $\mathrm{su}(m)$ model~\eqref{tJgen}-\eqref{tJV}
is supersymmetric if its coefficients are related by Eq.~\eqref{SUSY}. When this is the case this
general model reduces to~\eqref{tJSUSY}, which is equivalent to the~$\mathrm{su}(1|m)$ supersymmetric spin
chain~\eqref{hatH}. It should be stressed that the coefficient of the charge interaction
term~$n_in_j$ in the supersymmetric $t$-$J$ Hamiltonian~\eqref{tJSUSY} must depend on~$m$ as
specified in the latter equation, a fact that does not seem to have been previously noted in the
literature.

Note that Eqs.~\eqref{tJSUSY}--\eqref{hatH} are also valid for~$m=1$. In this case the
Hamiltonian~\eqref{tJSUSY} is that of a free fermion system (the terms proportional
to~$\mathbf{T}_i\cdot\mathbf{T}_j$ and $n_in_j$ vanish identically), as first noted by Göhmann and
Wadati~\cite{GW95}. This fact was recently exploited in Refs.~\cite{CFGRT16,CFGR17} to evaluate
the entanglement entropy of the ground state of (translation-invariant) $\mathrm{su}(1|1)$ spin chains of
HS type.

In the rest of this work we shall be mainly concerned with the supersymmetric $t$-$J$
model~\eqref{tJSUSY} with
\begin{equation}\label{tKY}
  t(x)=t(\pi/N)^2\sin^{-2}(\pi x/N)\,,
\end{equation}
which is the~$\mathrm{su}(m)$ version of the original KY model. In this case~\cite{CP78,FG05}
\begin{equation}
  \label{sinsum}
  \sum_{i\ne j}\sin^{-2}\bigl(\pi(i-j)/N\bigr)=\frac N3\,(N^2-1)\,,
\end{equation}
so that
\begin{equation}\label{t0HS}
  t_0=\frac{t\pi^2}{3N^2}(N^2-1)\,.
\end{equation}
Hence the Hamiltonian~$\hat H_0$ of the equivalent $\mathrm{su}(1|m)$ supersymmetric spin chain can be
written as
\begin{equation}
  \label{HtJHS}
  \hat H_0=\frac{2t\pi^2}{N^2}\bigg[ H_{\mathrm{HS}}^{(1|m)}-
  \frac 16(N^2-1){\mathcal F}\bigg]\,,
\end{equation}
where
\begin{equation}
  \label{HHS}
  H_{\mathrm{HS}}^{(1|m)}=\frac12\sum_{i<j}\frac{1-P_{ij}^{(1|m)}}{\sin^2(\pi(i-j)/N)}
\end{equation}
is the Hamiltonian of the $\mathrm{su}(1|m)$ Haldane--Shastry spin chain~\cite{Ha88,Sh88,Ha93,HH94}. In
fact, the model~\eqref{HtJHS} was introduced by Kawakami in the early 90's~\cite{Ka92}. On the
other hand, when~$t(x)$ is proportional to Eq.~\eqref{tde} the Hamiltonian~\eqref{tJSUSY} reduces
to the $\mathrm{su}(m)$ version of the original (nearest-neighbor) \hbox{$t$-$J$} model
\begin{equation}
  \fl
  H_0=t\sum_{i}{\mathcal P}\bigg[-\sum_{\sigma}(c^\dagger_{i\sigma}c_{i+1,\sigma}+c^\dagger_{i+1,\sigma}c_{i\sigma})
  +2\mathbf{T}_i\cdot\mathbf{T}_{i+1}-\left(1-\tfrac1m\right)n_in_{i+1}\bigg]{\mathcal P},
  \label{tJSUSYnn}
\end{equation}
where $N+1\equiv1$. The equivalent $\mathrm{su}(1|m)$ supersymmetric chain Hamiltonian is given by
\begin{equation}\label{US}
  \hat H_0=t\sum_i\big(1-P_{i,i+1}^{(1|m)}\big)-2t{\mathcal F}\,,
\end{equation}
where $P_{N,N+1}^{(1|m)}\equiv P_{1N}^{(1|m)}$, which is essentially the $\mathrm{su}(1|m)$
Uimin--Lai--Sutherland model~\cite{Ui70,La74,Su75}.

Next, inspired by Inozemtsev's elliptic spin chain~\cite{In90}, we introduce a one-parameter
family of supersymmetric $\mathrm{su}(m)$ $t$-$J$ models~\eqref{tJSUSY} which smoothly interpolate between
the $\mathrm{su}(m)$ KY model~\eqref{tJSUSY}-\eqref{tKY} and the (periodic) nearest-neighbors $\mathrm{su}(m)$
$t$-$J$ model~\eqref{tJSUSYnn}. More precisely, let
\begin{equation}\label{ellt}
  t(x)=t\bigg(\frac\alpha\pi\bigg)^{\!\!2}\!\sinh^2(\pi/\alpha)\bigg(\wp_N(x)
  -\frac{2\tilde\eta_1}{\alpha^2}\bigg)\,,
\end{equation}
where
\[
  \wp_N(x)\equiv\wp(x;N/2,\mathrm i\alpha/2),\quad \tilde\eta_1\equiv\zeta(1/2;1/2,\mathrm i N/(2\alpha))
\]
and~$\alpha>0$. In the latter formulas~$\wp(x;\omega_1,\omega_3)$ and $\zeta(x;\omega_1,\omega_3)$ denote
respectively the Weierstrass elliptic function with half-periods $\omega_1,\omega_3$ and its
corresponding zeta function, defined by
\begin{eqnarray*}
  \zeta(z;\omega_1,\omega_3)&=\frac1z+\sum_{l,n\in{\mathbb Z}^2-\{(0,0)\}}\bigg[\frac1{z-2l\omega_1-2n\omega_3}
                      +\frac{z}{(2l\omega_1+2n\omega_3)^2}\bigg]\,,\\
  \wp(z;\omega_1,\omega_3)&=-\zeta'(z;\omega_1,\omega_3)\,.
\end{eqnarray*}
It can be shown~\cite{In05,FG14JSTAT} that when~$1\le x\le N-1$
\[
  \lim_{\alpha\to0+}t(x)=t(\delta_{1,x}+\delta_{-1,x}),
\]
while the~$\alpha\to+\infty$ limit of Eq.~\eqref{ellt} is Eq.~\eqref{tKY}. Note that the constant
$t_0$ for the function~\eqref{ellt} is given by
\[
  t_0=\frac{2t}{\pi^2}\sinh^2(\pi/\alpha)(\tilde\eta_1-\eta_1)\,,\quad
  \eta_1\equiv\zeta(1/2;1/2,\mathrm i/(2\alpha))\
\]
(see, e.g., Ref.~\cite{FG14}). Thus the Hamiltonian~\eqref{hatH} of the~$\mathrm{su}(1|m)$ supersymmetric
spin chain equivalent to the~$\mathrm{su}(m)$ $t$-$J$ model~\eqref{tJSUSY} with elliptic
interactions~\eqref{ellt} is given by
\[
  \hat H_0 = tH_{I}^{(1|m)}-\frac{2t}{\pi^2}\,\sinh^2(\pi/\alpha)(\tilde\eta_1-\eta_1){\mathcal F},
\]
where
\[
  H_{I}^{(1|m)}=\bigg(\frac\alpha\pi\bigg)^{\!\!2}\!\sinh^2(\pi/\alpha)
  \sum_{i<j}\bigg(\wp_N(i-j)-\frac{2\tilde\eta_1}{\alpha^2}\bigg)(1-P_{ij}^{(1|m)})
\]
is the~$\mathrm{su}(1|m)$ version of Inozemtsev's elliptic spin chain~\cite{In90}. In fact, for $m=1$ the
partition function and thermodynamics of the latter chain were derived in Ref.~\cite{FG14JSTAT},
and the entanglement entropy of its ground state was analyzed in Ref.~\cite{CFGRT16}.

\section{Partition function}\label{sec.PF}

In this section we shall compute in closed form the partition function of the~$\mathrm{su}(m)$
KY model~\eqref{tJSUSY}-\eqref{tKY} by exploiting its equivalence with the
$\mathrm{su}(1|m)$ spin chain Hamiltonian of Haldane--Shastry type~\eqref{HtJHS}-\eqref{HHS}. As a matter
of fact, we shall consider the more general Hamiltonian
\begin{equation}
  \label{HKYmu}
  H=H_0-\frac12\sum_{\sigma=1}^{m-1}h_\sigma(n^\sigma-n^m)-\mu_c\sum_\sigma n^\sigma\equiv H_0+H_1,
\end{equation}
where $H_0$ is given by Eqs.~\eqref{tJSUSY}-\eqref{tKY} and
\[
  n^\sigma\equiv\sum_in_{i\sigma}
\]
denotes the total number of fermions of type~$\sigma$. The last term in~$H_1$ is the chemical
potential of the fermions (or, equivalently, of the total electric charge), while the first one
has a natural interpretation as arising from the interaction with an external $\mathrm{su}(m)$ magnetic
field with strengths $h_1,\dots,h_{m-1}$. Indeed, for $m=2$ the term $-(h_1/2)(n^1-n^2)$ equals
$-h_1S^z$, where~$S^z$ is the $z$ component of the total spin operator. This is indeed the
contribution to the energy arising from the interaction with the magnetic field $h_1\mathbf e_z$
of a charged fermion (with gyromagnetic ratio~$g=2$, and unit mass and electric charge in natural
units). More generally, for arbitrary $m\ge2$ we have
\[
  n^\sigma-n^m=\sum_i(n_{i\sigma}-n_{im}),
\]
where the operators $\mathrm i(n_{k\sigma}-n_{km})$ generate the standard $\mathrm{su}(m)$ Cartan subalgebra at the
$k$-th site. By Eq.~\eqref{Xissp}, the $\mathrm{su}(1|m)$ spin chain Hamiltonian~$\hat H$ equivalent to
$H$ is~$\hat H=\hat H_0+\hat H_1$, where
\[
  \hat H_1=-\frac12\sum_{\sigma=1}^{m-1}h_\sigma({\mathcal N}_\sigma-{\mathcal N}_m)-\mu_c{\mathcal F}
\]
and~${\mathcal N}_\sigma\equiv\sum_i\hat n^\sigma$ is the total numbers of fermions (created
by~$f^\dagger_{i\sigma}f_{i\sigma}$) of type~$\sigma$. More explicitly, we have
\begin{equation}
  \label{hatH0H1}
  \hat H=J H_{\mathrm{HS}}^{(1|m)}-\frac12\sum_{\sigma=1}^{m-1}h_\sigma({\mathcal N}_\sigma-{\mathcal N}_m)-
  (t_0+\mu_c){\mathcal F}\,,
\end{equation}
with $t_0$ given by Eq.~\eqref{t0HS} and $J=2t\pi^2/N^2$. The latter equation can be more
concisely rewritten as
\begin{equation}
  \label{hatHmu}
  \hat H=J H_{\mathrm{HS}}^{(1|m)}-\sum_\sigma\mu_\sigma{\mathcal N}_\sigma,
\end{equation}
where $\mu_\sigma$ is the chemical potential of the fermion of type~$\sigma$, given by
  \begin{eqnarray}\label{musi}
    \mu_\sigma&=\frac12\,h_\sigma+\mu_c+t_0\,,\quad 1\le\sigma\le m-1\,;\\
    \mu_m&=-\frac12\,\sum_{\sigma=1}^{m-1}h_\sigma+\mu_c+t_0\,.
           \label{musim}
  \end{eqnarray}

  From the above remarks, it follows that the partition function~${\mathcal Z}$ of the supersymmetric
  $\mathrm{su}(m)$ KY model~\eqref{HKYmu} coincides with that of the $\mathrm{su}(1|m)$ supersymmetric
  Haldane--Shastry chain with a chemical potential term given in Eq.~\eqref{hatHmu}. The partition
  function of the latter model has been recently evaluated in Ref.~\cite{FGLR18} by taking
  advantage of its connection with the supersymmetric spin Sutherland model via Polychronakos's
  freezing trick~\cite{Po93,Po94}, with the result
\begin{equation}\label{PF}
  {\mathcal Z}(q;{\boldsymbol{\mu}})=\sum_{\mathpalette\mathclapinternal{\mathbf{k}\in{\mathcal P}_N}}d(\mathbf{k})q^{\sum\limits_{i=1}^{r-1}JK_i(N-K_i)}
  \prod_{i=1}^{\mathpalette\mathclapinternal{N-r}}\big(1-q^{JK'_i(N-K'_i)}\big)\,.
\end{equation}
Here~$q\equiv\mathrm e^{-1/T}$, $\mathbf{k}=(k_1,\dots,k_r)\in({\mathbb N}\cup\{0\})^r$, ${\mathcal P}_N$ is the set of
partitions (with order taken into account) of the integer~$N$, $K_i\equiv\sum_{j=1}^ik_j$, and
the~$N-r$ positive integers~$K'_i$ are defined by
\[
  \{K_1',\dots,K'_{N-r}\}=\{1,\dots,N-1\}-\{K_1,\dots,K_{r-1}\}\,.
\]
For each multiindex~$\mathbf{k}\in({\mathbb N}\cup\{0\})^r$, the coefficient~$d(\mathbf{k})$ is defined by
\begin{equation}\label{dbk}
  d(\mathbf{k})=\prod_{i=1}^r\sum_{j=0}^{\min(k_i,m)}e_j(q^{-{\boldsymbol{\mu}}})\,,\quad
  q^{-{\boldsymbol{\mu}}}\equiv(q^{-\mu_1},\dots,q^{-\mu_m}),
\end{equation}
where
\begin{equation}\label{ejdef}
  e_j(y_1,\dots,y_m)\equiv\enspace\sum_{\mathpalette\mathclapinternal{1\le i_1<\cdots<i_j\le m}}\enspace y_{i_1}\cdots y_{i_j}
\end{equation}
denotes the elementary symmetric polynomial of degree~$j\le m$ in $m$ variables~$y_1,\dots,y_m$
(with~$e_0\equiv1$). By the above remarks, Eq.~\eqref{PF} yields also the partition function of
the supersymmetric~$\mathrm{su}(m)$ KY model~\eqref{HKYmu}. This is one of the main results of the present
work.

Equation~\eqref{dbk} for $d(\mathbf{k})$ can be considerably simplified introducing the
numbers~$\nu_i(\mathbf{k})$ defined by
\begin{eqnarray*}
  \nu_l(\mathbf{k})&=\left|\big\{i\in\{1,\dots,r\}:k_i=l\big\}\right|,\quad l=1,\dots,m-1\,,\\
  \nu_{m}(\mathbf{k})&=\left|\big\{i\in\{1,\dots,r\}:k_i\ge m\big\}\right|,
\end{eqnarray*}
where~$|A|$ denotes the cardinal of the set~$A$. We then have
\begin{equation}
  \label{degm}
  d(\mathbf{k})=\prod_{l=1}^{m-1}\Bigg[\sum_{j=0}^le_j(q^{-{\boldsymbol{\mu}}})\Bigg]^{\nu_l(\mathbf{k})}\cdot
  \prod_{j=1}^m(1+q^{-\mu_j})^{\nu_m(\mathbf{k})}.
\end{equation}
For instance, in the $\mathrm{su}(2)$ and~$\mathrm{su}(3)$ cases we respectively have
\[
  d(\mathbf{k})=(1+q^{-\mu_1}+q^{-\mu_2})^{\nu_1(\mathbf{k})}[(1+q^{-\mu_1})(1+q^{-\mu_2})]^{\nu_2(\mathbf{k})}
\]
and
\[
  \fl
  d(\mathbf{k})=\Big(1+{\textstyle\sum\limits_\sigma}\,q^{-\mu_\sigma}\Big)^{\nu_1(\mathbf{k})}
  \Big(1+{\textstyle\sum\limits_\sigma}\,q^{-\mu_\sigma}
  +{\textstyle\sum\limits_{\sigma<\sigma'}}q^{-(\mu_\sigma+\mu_{\sigma'})}\Big)^{\nu_2(\mathbf{k})}
  \prod_\sigma(1+q^{-\mu_\sigma})^{\nu_3(\mathbf{k})}\,.
\]
The latter expressions becomes even simpler when applied to the ``pure'' supersymmetric KY model
(i.e., without magnetic field or chemical potential terms), for which~$\mu_\sigma=t_0$ for
all~$\sigma=1,\dots,m$ according to Eqs.~\eqref{musi}-\eqref{musim}. Indeed, in this case we have
\[
  e_j(q^{-{\boldsymbol{\mu}}})=e_j(q^{-t_0},\dots,q^{-t_0})=\sum_{\mathpalette\mathclapinternal{1\le i_1<\cdots<i_j\le m}}\enspace
  q^{-jt_0}=\binom mj q^{-jt_0}\,,
\]
and hence
\[
  d(\mathbf{k})=\big(1+q^{-t_0}\big)^{m\nu_m(\mathbf{k})}\prod_{l=1}^{m-1}\Bigg[\sum_{j=0}^l
  \binom mj q^{-jt_0}\Bigg]^{\nu_l(\mathbf{k})}.
\]
For instance, the partition function of the original ($\mathrm{su}(2)$) supersymmetric KY model is given
by
\[
  \fl
  {\mathcal Z}_N(q)=\sum_{\mathpalette\mathclapinternal{\mathbf{k}\in{\mathcal P}_N}}\big(1+2q^{-t_0}\big)^{\nu_1(\mathbf{k})}\big(1+q^{-t_0}\big)^{2\nu_2(\mathbf{k})}
  q^{\sum\limits_{i=1}^{r-1}JK_i(N-K_i)}
  \prod_{i=1}^{\mathpalette\mathclapinternal{N-r}}\big(1-q^{JK'_i(N-K'_i)}\big)\,.
\]

\section{Spectrum and motifs}\label{sec.motifs}

In this section we shall give a complete description of the spectrum of the supersymmetric
$\mathrm{su}(m)$ KY model~\eqref{tJSUSY}-\eqref{tKY} ---or, more generally, the
Hamiltonian~\eqref{HKYmu}--- in each subspace with well-defined spin content in terms of the
supersymmetric version of Haldane's motifs~\cite{Ha93} and their associated skew Young
tableaux~\cite{BBHS07,BBH10,FG15}. In particular, this description implies the validity of the
Saiga--Kuramoto conjecture, which is one of the main results of this paper.

We start by recalling that the partition function~\eqref{PF} of the~$\mathrm{su}(1|m)$ spin
chain~\eqref{hatH0H1} ---and, hence, of the~supersymmetric $\mathrm{su}(m)$ KY model~\eqref{HKYmu}---
exactly coincides with the partition function of the inhomogeneous vertex model with
energies~\cite{BBH10,FGLR18}
\begin{equation}\label{specH}
  E(\mathbf{s})=J\sum_{i=1}^{N-1}\delta(s_i,s_{i+1})i(N-i)
  -\sum_{i}\mu_{s_i}\,,
\end{equation}
where~$\mu_0\equiv0$, $\mathbf{s}\in\{0,\dots,m\}^N$ and~$\delta(s,s')$ is defined by
\[
  \delta(s,s')=\cases{1,&$s>s' \enspace\text{or}\enspace s=s'>0$\\
    0,&$s<s'\enspace\text{or}\enspace s=s'=0$\,.
  }
\]
The first sum in Eq.~\eqref{specH} can be interpreted as the energy of a one-dimensional vertex
model with~$N+1$ vertices~$0,\dots,N$ joined by~$N$ bonds with
values~$s_1,\dots,s_N\in\{0,\dots,m\}$, the energy associated to the $i$-th vertex~being equal
to~$\delta(s_i,s_{i+1})i(N-i)$. For this reason, we shall henceforth refer to the vector~$\mathbf{s}$ as
the \emph{bond vector}. Likewise, the vectors~$\boldsymbol{\delta}(\mathbf{s})$ with components~$\delta(s_i,s_{i+1})$
$(1\le i\le N-1)$ in Eq.~\eqref{specH} can be identified with $\mathrm{su}(1|m)$
motifs~\cite{Ha93,HB00,BBH10}. Thus the spectrum of the $\mathrm{su}(m)$ KY model (with the correct
degeneracy for each level) can be computed from Eq.~\eqref{specH} by letting~$\mathbf{s}$ run over all
possible~$(m+1)^N$ bond vectors. It is important to note that the energies~\eqref{specH} depend not
only on the motif~$\boldsymbol{\delta}$ but also on the chemical potentials~$\mu_\alpha$ through the last term. This
term will in general break the huge degeneracy associated to the motifs~$\boldsymbol{\delta}$, which is in part
due to the invariance of the model~\eqref{hatHmu} with $\mu_\alpha=0$ (i.e., the $\mathrm{su}(1|m)$
supersymmetric HS chain) under the Yangian~$Y(\mathrm{gl}(1|m))$~\cite{FG15}. In other words, the general
model~\eqref{hatHmu} should be far less degenerate than the $\mathrm{su}(1|m)$ supersymmetric HS spin
chain.

Let us denote by $\psi_\mathbf{s}$ the unique eigenfunction of the~supersymmetric $\mathrm{su}(m)$ KY
model~\eqref{HKYmu} corresponding to the eigenvalue $E(\mathbf{s})$ \eqref{specH} associated with the
bond vector~$\mathbf{s}$. Our aim is to determine the magnon numbers or spin content of the
eigenfunction $\psi_\mathbf{s}$ directly from the structure of its bond vector~$\mathbf{s}$. To this end, we
shall now show how the spectrum of the restriction of the Hamiltonian~\eqref{HKYmu} to each
subspace with well-defined spin content can be fully generated from the motif
formula~\eqref{specH} by suitably restricting the components of the bond vector~$\mathbf{s}$.

More precisely, let us denote by
\[
  n(N_0,\dots,N_m)\equiv n(\mathbf N)\,,\qquad N_0+\cdots+N_m\equiv|\mathbf N|=N,
\]
the subspace of the Hilbert space of the~$\mathrm{su}(m)$ KY model in which each number operator~$n^{\alpha}$
has a well-defined value $N_\alpha$. We shall henceforth refer to the vector~$\mathbf N$ as the
\emph{magnon content} of the subspace~$n(\mathbf N)$. Note next that, by Eqs.~\eqref{HKYmu}
and~\eqref{hatH0H1}, the KY Hamiltonian~$H$ in Eq.~\eqref{HKYmu} can be written as
\begin{equation}
  H= {H}_0'-\sum_\sigma \mu_\sigma n^\sigma \,,\qquad
  H_0' \equiv H_0+t_0\sum_\sigma n^\sigma\,,
 \label{HKYmu1}
\end{equation}
where $ H_0$ and $t_0$ are respectively given by Eqs.~\eqref{tJSUSY}-\eqref{tKY} and~\eqref{t0HS}.
Clearly, $H_0'$ does not depend on the chemical potentials~$\mu_\sigma$ (since~$H_0$ is also
independent of~$\mu_\sigma$); in fact, $\widehat {H_0'} =J H_{\mathrm{HS}}^{(1|m)}$.
It should be noted at this point that the subspaces $n(\mathbf N)$ are invariant under both $H$ and
${H}_0'$, since the number operators $n^\sigma$ obviously commute with the
operators~$c^\dagger_{k\sigma}c_{l\sigma}$, $n_k$, $\mathbf{T}_k\cdot\mathbf{T}_l$.
Consequently, by Eq.~\eqref{HKYmu1}, the partition function of the Hamiltonian~$H$ can be written
as
\begin{equation}\label{Z0bN}
  {\mathcal Z}(q;{\boldsymbol{\mu}})=\sum_{\mathbf N;|\mathbf N|=N}q^{-\sum_\sigma \mu_\sigma N_\sigma}{\mathcal Z}_0^{\mathbf N}(q),
\end{equation}
where~${\mathcal Z}_0^{\mathbf N}(q)$ denotes the partition function of the restriction of~${H}_0'$ to the
subspace~$n(\mathbf N)$. We also know that
\begin{equation}\label{equiv}
  {\mathcal Z}(q;{\boldsymbol{\mu}})={\mathcal Z}_V(q;{\boldsymbol{\mu}})\,,
\end{equation}
where~${\mathcal Z}_V(q;{\boldsymbol{\mu}})$ denotes the partition function of the vertex model~\eqref{specH}. Let us now
rewrite the latter equation as
\[
  E(\mathbf{s})=J\sum_{i=1}^{N-1}\delta(s_i,s_{i+1})i(N-i) -\sum_{\sigma}\mu_\sigma N_{\sigma}(\mathbf{s}),
\]
where~$N_\sigma(\mathbf{s})$ denotes the number of components of the bond vector~$\mathbf{s}$ equal to~$\sigma$.
We then have
\begin{equation}\label{ZNsi}
  \fl
  {\mathcal Z}_V(q;{\boldsymbol{\mu}})=\sum_{\mathbf{s}}q^{E(\mathbf{s})}=\sum_{\mathpalette\mathclapinternal{\mathbf N;|\mathbf N|=N}}q^{-\sum_\sigma\mu_\sigma
    N_\sigma}\sum_{\mathpalette\mathclapinternal{\mathbf{s};N_\alpha(\mathbf{s})=N_\alpha}}q^{J\sum_{i=1}^{N-1}\delta(s_i,s_{i+1})i(N-i)},
\end{equation}
where the last sum is restricted to all bond vectors~$\mathbf{s}$ such that~$N_\alpha(\mathbf{s})=N_\alpha$ for
all~$\alpha=0,\dots,m$. Note that the right-hand sides of Eqs.~\eqref{Z0bN} and~\eqref{ZNsi} are both
polynomials in the variables~$x_\alpha\equiv\mathrm e^{\beta\mu_\alpha}$ ($\alpha=0,\dots,m$).
Equating the coefficient of~$x_0^{N_0}\cdots x_m^{N_m}$ in both equivalent expressions
for~${\mathcal Z}(q;{\boldsymbol{\mu}})$ we immediately deduce that the partition function of the restriction of $H_0'$
to the subspace~$n(\mathbf N)$ is given by
\begin{equation}
  \label{Z0motifs}
  {\mathcal Z}_0^{\mathbf N}(q)=\sum_{\mathbf{s}:N_\alpha(\mathbf{s})=N_\alpha}q^{J\sum_{i=1}^{N-1}\delta(s_i,s_{i+1})i(N-i)}\,.
\end{equation}
It follows from the latter equation that the spectrum of~$H_0'$ in the subspace~$n(\mathbf N)$ can be
generated from the motif formula
\[
  \fl E_0(\mathbf{s};\mathbf N)=J\sum_{i=1}^{N-1}\delta(s_i,s_{i+1})i(N-i),\qquad \text{with}\enspace
  N_\alpha(\mathbf{s})=N_\alpha\enspace\text{for all }\alpha=0,\dots,m\,.
\]
In view of Eq.~\eqref{HKYmu1}, the spectrum of the full Hamiltonian~$H$ restricted to the
subspace~$n(\mathbf N)$ (with the correct degeneracies for all levels) is generated by the analogous
formula
\begin{equation}
  \label{specHcNmotifs}
  \fl
  E(\mathbf{s};\mathbf N)=E_0(\mathbf{s};\mathbf N)-\sum_\sigma\mu_\sigma N_\sigma,\qquad
  \text{with}\enspace N_\alpha(\mathbf{s})=N_\alpha\enspace\text{for all }\alpha=0,\dots,m\,.
\end{equation}
We thus conclude that the spectrum of~$H$ on~$n(\mathbf N)$ is obtained from Eq.~\eqref{specH} by
imposing the natural conditions~$N_\alpha(\mathbf{s})=N_\alpha$ on the bond vector~$\mathbf{s}$. It follows from
this assertion that we can label the eigenfunctions of~$H$ in such a way that each $\psi_\mathbf{s}$
belongs to the subspace~$n(\mathbf N(\mathbf{s}))$ containing exactly~$N_\alpha(\mathbf{s})$ magnons of type~$\alpha$.
This is in fact one of the main results of this section, whose consequences we shall explore next.

To begin with, we shall use the method of Ref.~\cite{BBH10} to express the
spectrum~\eqref{specHcNmotifs} of the supersymmetric KY model~\eqref{HKYmu} on the invariant
subspace $n(\mathbf N)$ in an alternative way. To this end, we first note that the numerical values of
the energies can be computed from the formula
\begin{equation}\label{EdeN}
  E_{\boldsymbol{\delta},\mathbf N}=J\sum_{i=1}^{N-1}\delta_i\cdot i(N-i)-\sum_\sigma\mu_\sigma N_\sigma\,,
\end{equation}
where now~$\boldsymbol{\delta}\equiv(\delta_1,\dots,\delta_{N-1})$ is a supersymmetric motif, i.e., a sequence of~$N-1$
zeros and ones. Secondly, the degeneracy of each energy~$E_{\boldsymbol{\delta},\mathbf N}$ (which could be zero) is
evaluated by counting the number of ways of filling the border strip associated to the
motif~$\boldsymbol{\delta}$ according to the usual~$\mathrm{su}(1|m)$ rules~\cite{KKN97,FG15}, with the additional
restriction that each number~$\alpha\in\{0,\dots,m\}$ must appear exactly~$N_\alpha$ times. More
precisely, given the motif $(\delta_1,\dots,\delta_{N-1})$ its associated border strip is constructed by
starting with one box, and then reading the motif from left to right and adding a box below
(resp.~to the left of) the $i$-th box provided that $\delta_i$ is equal to $0$ (resp.~$1$); see,
e.g., Fig.~\ref{fig.yd}. This border strip should then be filled with the numbers $0,1,\dots,m$
according to the following rules:
\begin{enumerate}[i)]
\item The numbers in each row form a nondecreasing sequence, allowing only the repetition of
  positive numbers.
\item The numbers in each column (read from top to bottom) form a nondecreasing sequence,
  allowing only the repetition of $0$.
\item Each number~$\alpha\in\{0,\dots,m\}$ must appear exactly~$N_\alpha$ times.
\end{enumerate}
Each filling of a border strip according to the previous rules is called a (skew) Young tableau.
Given such a tableau, it is straightforward to check that its associated motif can be obtained
from the bond vector~$(s_1,\dots,s_N)$ whose components are the numbers in the tableau read from
top to bottom by setting~$\delta_i=\delta(s_i,s_{i+1})$. The equivalence between the description of the
spectrum through Eq.~\eqref{EdeN} (where the degeneracy of a motif~$\boldsymbol{\delta}$ is evaluated counting
all the fillings of its associated border strip allowed by rules i)--iii)) and
Eq.~\eqref{specHcNmotifs} essentially follows from this observation. Note also that, by
Eq.~\eqref{EdeN}, all the Young tableaux associated to a given motif have the same energy within
each invariant subspace~$n(\mathbf N)$ .
\begin{figure}[t]
  \hfill\small
\[
  \ytableaushort{\none0,01,0,2
}
\qquad
  \ytableaushort{\none0,02,0,1
}
\qquad
  \ytableaushort{\none1,02,0,0
}
\]
\hfill
\caption{Allowed Young tableaux for the motif~$(0,1,0,0)$ in the subspace~$n(3,1,1)$.
    \label{fig.yd}}
\end{figure}

Consider, for instance, the motif~$\boldsymbol{\delta}=(0,1,0,0)$ for $N=5$ particles in the~$\mathrm{su}(1|2)$ case. As
explained above, the degeneracy associated to this motif in each invariant subspace~$n(\mathbf N)$ is
given by all possible Young tableaux generated from it according to the above three rules. For
instance, for the subspace~$n(3,1,1)$ it is easy to check that there are exactly three allowed
Young tableaux for the above motif (cf.~Fig.~\ref{fig.yd}), with energy $6J-\mu_1-\mu_2$. In fact,
it is a straightforward matter to verify that there are exactly~$8$ invariant subspaces~$n(\mathbf N)$
with fixed spin content compatible with the motif~$(0,1,0,0)$ (i.e., with nonzero degeneracy),
whose respective degeneracies and energies are listed in Table~\ref{table.subs}.
\begin{table}[h]
  \centering
  \begin{tabular}{|c|c|c|}\hline
    $\mathbf N$& Degeneracy& Energy\\
    \hline
    (1,2,2)& 1& $6J-2(\mu_1+\mu_2)$\\
    (3,0,2)& 1& $6J-2\mu_2$\\
    (3,2,0)& 1& $6J-2\mu_1$\\
    (4,0,1)& 1& $6J-\mu_2$\\
    (4,1,0)& 1& $6J-\mu_1$\\
    (2,1,2)& 2& $6J-(\mu_1+2\mu_2)$\\
    (2,2,1)& 2& $6J-(2\mu_1+\mu_2)$\\
    (3,1,1)& 3& $6J-(\mu_1+\mu_2)$\\
    \hline
  \end{tabular}
  \caption{Invariant subspaces~$n(\mathbf N)$ compatible with the motif~$(0,1,0,0)$ in the $\mathrm{su}(1|2)$
    case.}
  \label{table.subs}
\end{table}

The above general results considerably simplify for the supersymmetric $\mathrm{su}(m)$ KY model without
magnetic field or chemical potential terms, given by Eq.~\eqref{tJSUSY}-\eqref{tKY}. Indeed, in
this case $h_\sigma=\mu_c=0$, so that from Eq.~\eqref{musim} we obtain
\[
  \mu_{\sigma}=t_0=\frac{t\pi^2}{3N^2}(N^2-1),\qquad 1\le\sigma\le m\,.
\]
Thus Eq.~\eqref{EdeN} becomes
\begin{equation}\label{EKY}
  \fl
  E_{\boldsymbol{\delta},\mathbf N}=J\sum_{i=1}^{N-1}\delta_i\cdot i(N-i)-t_0(N-N_0)
  =J\left[\sum_{i=1}^{N-1}\delta_i\cdot i(N-i)-\frac16\,(N^2-1)(N-N_0)\right],
\end{equation}
which depends on~the spin content~$\mathbf N$ only through~$N_0$. Thus for a given motif~$\boldsymbol{\delta}$ all of
its compatible invariant subspaces~$n(\mathbf N)$ with the same number of holes~$N_0$ will now have the
same energy. For instance, for the case~$N=5$, $m=2$ and the motif~$(0,1,0,0)$ considered above
there are only four sectors with different energies, corresponding to $N_0=1$ (singlet), $2$ (four
times degenerate), $3$ (five times degenerate) and~$4$ (twice degenerate), with respective
energies~$-2J(7-2N_0)$. In general, in order to compute the degeneracy associated to each
motif~$\boldsymbol{\delta}$ in a sector with a given number of holes~$N_0$ we just have to count the number of
allowed Young tableaux according to rules i) and~ii) above, replacing rule three by the following:

\begin{description}
\item{iii')\enspace} The number~$0$ must appear exactly~$N_0$ times.
\end{description}
The rules i), ii), iii') provide a complete description of the spectrum of the
supersymmetric~$\mathrm{su}(m)$ KY model~\eqref{tJSUSY}-\eqref{tKY}, with the correct degeneracy for each
energy~\eqref{EKY}, for arbitrary values of $m$ and the number of sites~$N$. In particular, when
applied to the simplest $\mathrm{su}(2)$ case these rules provide the first rigorous proof known to the
authors to the long-standing conjecture of Saiga and Kuramoto mentioned in the Introduction.
\begin{table}[t] \centering
  \begin{tabular}{|c|c|c|c|c|c|c|c|}\hline
    Motif &0 &1 &2 &3 &4 &5 & Energy\\
    \hline
    $(0,0,0,0)$ &0 &0 &0 &1 &2 &1 & $-4J(5-N_0)$\\
    $(0,0,0,1$), $(1,0,0,0)$ &0 &0 &2 &4 &2 &0& $-4J(4-N_0)$\\
    $(0,0,1,0)$, $(0,1,0,0)$ &0 &1 &4 &5 &2 &0 & $-2J(7-2N_0)$\\
    $(1,0,0,1)$ &0 &4 &8 &4 &0 &0 &$-4J(3-N_0)$\\
    $(0,0,1,1)$, $(1,1,0,0)$ &0 &3 &6 &3 &0 &0 & $-2J(5-2N_0)$\\
    $(1,0,1,0)$, $(0,1,0,1)$ &2 &8 &10 &4 &0 &0 & $-2J(5-2N_0)$\\
    $(0,1,1,0)$ &2 &7 &8 &3 &0 &0 & $-4J(2-N_0)$\\
    $(1,0,1,1)$, $(1,1,0,1)$ &6 &12 &6 &0 &0 &0 & $-2J(3-2N_0)$\\
    $(0,1,1,1)$, $(1,1,1,0)$ &4 &8 &4 &0 &0 &0 & $-4J(1-N_0)$\\
    $(1,1,1,1)$ &6 &5 &0 &0 &0 &0 & $4JN_0$\\
    \hline
  \end{tabular}
  \caption{Spectrum of the $\mathrm{su}(1|2)$ KY model~\eqref{tJSUSY}-\eqref{tKY} with~$N=5$ sites. The
    integers appearing in the columns labeled $0, \dots, 5$ are the degeneracies of the motif(s)
    in each row for the subspaces with $N_0=0, \dots, 5$ holes.}
  \label{table.spec}
\end{table}

As a simple illustration of the above assertion, we present in Table~\ref{table.spec} the detailed
spectrum of the (original) $\mathrm{su}(2)$ KY model~\eqref{tJSUSY}-\eqref{tKY} with~$N=5$ sites. More
precisely, in each of the columns of this table labeled with the integers $N_0=0,\dots,5$ we list
the degeneracy associated to the motif(s) in each row for the subspace with $N_0$ holes. This
degeneracy is computed by first generating all the skew Young tableaux compatible with each of the
$2^{N-1}=16$ possible motifs according to rules i), ii) and iii') above, which can be easily
accomplished using a simple \emph{Mathematica} program. We then sort the resulting tableaux
according to the number of holes (zeros) in each of them. By Eq.~\eqref{EKY}, a
motif~$\boldsymbol{\delta}=(\delta_1,\dots,\delta_{N-1})$ and its reverse~$\boldsymbol{\delta}'=(\delta_{N-1},\dots,\delta_1)$ clearly have
the same energy. In fact, it can be shown that two such motifs give rise to the same numbers of
compatible Young tableaux in each subspace~$n(\mathbf N)$ (see the appendix for details). For this
reason, we have grouped together in Table~\ref{table.spec} two motifs that are the reverse of each
other. More generally, if we exchange any two components~$\delta_k$ and~$\delta_{N-k}$ of a motif~$\boldsymbol{\delta}$
we obtain a motif with the same energy as~$\boldsymbol{\delta}$ in each subspace with~$N_0$ holes. However, these
two motifs may not necessarily yield the same number of compatible Young tableaux in such a
subspace. For instance, the motifs $(0,0,1,1)$ and $(1,0,1,0)$ have the same energy~$-2J(5-2N_0)$
in a subspace with~$N_0$ holes, but it is clear from Table~\ref{table.spec} that their
degeneracies differ in each of these subspaces for $N_0=0,\dots,3$. Finally, it is apparent from
Table~\ref{table.spec} that the ground state is obtained from the motifs~$(0,0,1,0)$, $(0,1,0,0)$
(in the sector with one hole) and $(1,0,1,0)$, $(0,1,0,1)$ (in the sector with no holes). It is
thus six times degenerate, with energy $-10J$.

\section{Ground state phases for the supersymmetric spin~\mbox{\boldmath$1/2$} KY model}
\label{sec.GS}

The complete description of the spectrum of the supersymmetric KY model in terms of motifs, bond
vectors and their associated skew Young tableaux developed in the previous section is particularly
suited to studying its ground state. As an example, we shall compute next the ground state energy
per site of the~spin~$1/2$ KY model in the thermodynamic limit for all possible values of the
magnetic field strength~$h\equiv h_1$ and chemical potential~$\mu\equiv\mu_c+t_0$. To this end,
let us first choose the unit of energy so that $t=1/(2\pi^2)$, so that Eq.~\eqref{EdeN} reads
\begin{equation}\label{Ensu2}
E_{\boldsymbol{\delta},\mathbf N}=\frac1{N^2}\sum_{i=1}^{N-1}\delta_i\cdot i(N-i)-\frac h2\,(N_1-N_2)-\mu(N_1+N_2)\,.
\end{equation}
We shall also assume without loss of generality that~$h\ge0$, since taking $h<0$ simply reverses
the role of the ``up'' ($\sigma=1$) and ``down'' ($\sigma=2$) fermions. It is then clear from the term
proportional to~$h$ in Eq.~\eqref{Ensu2} that the ground state(s) must belong to an invariant
subspace~$n(\mathbf N)$ with~$N_1\ge N_2$. Since the dispersion function~$i(N-i)$ is symmetric
about~$i=N/2$ and has an absolute maximum at this point, by Eq.~\eqref{specH} for $N$ large enough
the bond vector minimizing the energy in the subspace~$n(\mathbf N)$ must be of the form\footnote{We are
  actually assuming here that both $N_1$ and~$N_2$ are even. In other cases the form of the
  minimizing bond vector differs slightly from Eq.~\eqref{de0}, but the formula for its energy
  coincides in the thermodynamic limit with the one given below. Note also that unless $N_1$
  and~$N_2$ are both odd there is actually an additional bond vector with the same energy
  as~\eqref{de0} (or its variants, for~$N_1$ and $N_2$ of opposite parity).}
\begin{equation}\label{de0}
  \mathbf{s}_0=(\,\underbrace{1\cdots1}_{\mathpalette\mathclapinternal{(N_1-N_2)/2}}\,\underbrace{12\cdots12}_{N_2}\,
  \underbrace{0\cdots0}_{N_0}\,\underbrace{12\cdots12}_{N_2}\,
  \underbrace{1\cdots1}_{\mathpalette\mathclapinternal{(N_1-N_2)/2}}\,)\,.
\end{equation}
By Eq.~\eqref{Ensu2}, the energy per site of the corresponding motif
\[
   \boldsymbol{\delta}_0=(\,\underbrace{1,\dots,1,}_{\mathpalette\mathclapinternal{(N_1-N_2)/2}}\,\underbrace{0,1,\dots,0,1,}_{N_2}\,
   \underbrace{0,\dots,0,}_{N_0}\,\underbrace{0,1,\dots,0,1,}_{N_2}\,
   \underbrace{1,\dots,1}_{\mathpalette\mathclapinternal{(N_1-N_2)/2-1}}\,)
\]
is given by
\[
  \frac{E_{\boldsymbol{\delta}_0,\mathbf N}}N=\frac2N\sum_{k=1}^{Nt}\varepsilon(x_k)
  +\frac2N\sum_{k=1}^{N(s-t)/2}\kern-6pt\varepsilon(x_{Nt+2k})
  -ht-2\mu s\equiv u_{\boldsymbol{\delta}_0,\mathbf N}\,,
\]
where the factor of $2$ before the sums comes from the obvious symmetry of~$\boldsymbol{\delta}_0$ around $N/2$
and we have set
\[
  s=\frac{N_1+N_2}{2N},\qquad t=\frac{N_1-N_2}{2N},\qquad x_k=\frac{k}{N}\,,
  \qquad \varepsilon(x)=x(1-x)\,.
\]
Hence in the thermodynamic limit the minimum energy per site in a subspace~$n(\mathbf N)$ with
$N_1\ge N_2$ is given by
\begin{eqnarray}
  \fl
  u(s,t)\equiv\lim_{N\to\infty}u_{\boldsymbol{\delta}_0,\mathbf N}
  &=2\int_0^t\varepsilon(x)\,\mathrm{d} x+\int_t^s\varepsilon(x)\,\mathrm{d}
    x-ht-2\mu s\nonumber\\
  \fl
  &=\int_0^s\varepsilon(x)\,\mathrm{d} x+\int_0^t\varepsilon(x)\,\mathrm{d} x-ht-2\mu s
    \equiv f_{2\mu}(s)+f_h(t)\,,
    \label{ude0}
\end{eqnarray}
where
\[
  f_\lambda(s)=\int_0^s\varepsilon(x)\,\mathrm{d} x-\lambda s\equiv f(s)-\lambda s\,.
\]
The ground state energy of the spin~$1/2$ KY model in the thermodynamic limit is thus the minimum
value of the function $u(s,t)$ in the triangle
\[
  D=\{(s,t)\in{\mathbb R}^2:0\le t\le s\le 1/2\}\,.
\]
In order to compute this minimum value, note first of all that~$f'_\lambda(x)=\varepsilon(x)-\lambda$, with
$\varepsilon(x)$ monotonically increasing from~$0$ to $1/4$ in the interval~$[0,1/2]$. Hence $f_\lambda$ is
monotonically increasing\footnote{For the sake of conciseness, we shall implicitly assume in what
  follows that the functions~$f_\lambda$ and~$\varepsilon$ are restricted to the interval of
  interest~$[0,1/2]$.} over the interval $[0,1/2]$ for~$\lambda\le0$, monotonically decreasing
for~$\lambda\ge1/4$, and has a unique global minimum at the point
\[
  x_0(\lambda)\equiv\varepsilon^{-1}(\lambda)=\frac12\,\big(1-\sqrt{1-4\lambda}\,\big)\in(0,1/2)
\]
for $0<\lambda<1/4$. We thus have the following possibilities (recall that we are assuming throughout
that $h\ge0$):

\medskip

\noindent i)\enspace $\displaystyle\frac h2<\mu<1/8$

\smallskip
\noindent In this case~$f_h(t)$ and $f_{2\mu}(s)$ both have a unique minimum over the
interval~$[0,1/2]$ respectively at the points~$t_0=x_0(h)\in[0,1/2)$
and~$s_0=x_0(2\mu)\in(0,1/2)$, with~$t_0< s_0$ since~$\varepsilon(t_0)=h<2\mu=\varepsilon(s_0)$. Hence~$u$
attains its global minimum on~$D$ at the point~$(s_0,t_0)$ (which is an interior point if~$h>0$),
and the ground state energy density~$u$ is consequently given by
\begin{eqnarray}
  u&=u(s_0,t_0)=\int_0^{x_0(h)}\varepsilon(x)\,\mathrm{d} x+\int_0^{x_0(2\mu)}\varepsilon(x)\,\mathrm{d} x-hx_0(h)-2\mu
     x_0(2\mu)\nonumber\\
   &=\frac1{6}\,\bigg[1-3(h+2\mu)-\frac12(1-4h)^{3/2}-\frac12(1-8\mu)^{3/2}\bigg]\,.
     \label{uT}
\end{eqnarray}
Note also that $t=(N_1-N_2)/(2N)$ and~$2s=(N_1+N_2)/N$ are respectively equal to the magnetization
and the charge density per site (assuming that the fermions have unit charge and gyromagnetic
ratio equal to~2). Hence the zero-temperature magnetization and charge densities are given by
\[
  m_s=x_0(h)=\frac12\,\big(1-\sqrt{1-4h}\,\big)\,,\qquad
  n_c=2x_0(2\mu)=1-\sqrt{1-8\mu}\,.
\]
The corresponding magnetic and charge susceptibilities are obtained by differentiation, namely,
\[
  \fl
  \chi_s=\pdf{m_s}h=(1-4h)^{-1/2}=(1-2m_s)^{-1}\,,\quad
  \chi_c=\pdf{n_c}\mu=4(1-8\mu)^{-1/2}=4(1-n_c)^{-1}\,,
\]
in agreement with known results (see, e.g.,~\cite{KK95}).

\medskip\noindent

If $(h,\mu)$ lies outside the region~$\{(h,\mu)\in{\mathbb R}^2:0\le h/2<\mu<1/8\}$, the system
\[
  \varepsilon(s)=2\mu\,,\qquad \varepsilon(t)=h
\]
determining the critical points of $u$ has no solutions within the interior of~$D$. Thus the
function~$u(s,t)$ must necessarily attain its minimum value in the triangle~$D$ on its sides. It
is worth noting in this respect that this minimum cannot be reached at an interior point of the
horizontal side~$t=0$ unless~$h=0$, since for~$h>0$ we have~$f_h'(0)=-h<0$ and thus
\[
  u(s,0)=f_{2\mu}(s)>f_{2\mu}(s)+f_h(t)=u(s,t)
\]
for sufficiently small~$t>0$. We are left with the following possibilities:

\medskip
\noindent ii)\enspace $h\ge1/4$, $\mu+\displaystyle\frac h2\ge1/4$

\smallskip\noindent In this case $u(s,s)=f_h(s)+f_{2\mu}(s)=2f_{h/2+\mu}(s)$
and~$u(1/2,t)=f_h(t)+f_{2\mu}(1/2)$ are both decreasing, so that these functions have a unique
global minimum on the interval~$[0,1/2]$ at the right endpoint~$1/2$, with common value
\begin{equation}\label{uminW1}
  u(1/2,1/2)=2f_{h/2+\mu}(1/2)=
  \frac1{6}\,\big[1-3(h+2\mu)\big]\,.
\end{equation}
This must be the unique minimum of~$u$ on~$D$, since as~$h>0$ this minimum cannot be attained on
the interior of the side~$t=0$. Hence in this case the ground state energy per site is given by
Eq.~\eqref{uminW1}, while the zero-temperature magnetization and charge densities are simply
\[
  n_c=2m_s=1\,.
\]
Note that~$s=t=1/2$ is equivalent to~$N_1=N$, so that this phase consists only of up fermions.

\medskip
\noindent iii)\enspace $\mu+\displaystyle\frac h2\le0$

\smallskip\noindent Since~$\mu\le -h/2\le0$, the minimum of~$u(s,0)=f_{2\mu}(s)$
and~$u(s,s)=2f_{h/2+\mu}(s)$ on~$[0,1/2]$ is located at~$s=0$, while
\[
  u(1/2,t)=f(t)+f(1/2)-ht-\mu> -\frac h2-\mu\ge 0=u(0,0)\,.
\]
Thus in this case the ground state energy, zero-temperature magnetization and charge all vanish.
In fact, since $s=t=0$ is equivalent to $N_0=N$, this is the trivial phase consisting only of
holes.

\medskip
\noindent iv)\enspace $h<1/4$, $\mu\ge1/8$

\smallskip\noindent In this case $u(1/2,t)$ attains its minimum value in the interval~$[0,1/2]$ at
$t_0=x_0(h)\in[0,1/2)$. On the other hand,
\[
  u(s,s)=f_h(s)+f_{2\mu}(s)> f_h(t_0)+f_{2\mu}(1/2)=u(1/2,t_0),
\]
since $f_{2\mu}$ is decreasing in the interval~$[0,1/2]$ on account of the condition~$\mu\ge1/8$.
Similarly,
\[
  u(s,0)=f_{2\mu}(s)\ge f_{2\mu}(1/2)=u(1/2,0)\ge u(1/2,t_0),
\]
with equality only if~$h=0$ and~$s=1/2$. Thus the global minimum of~$u$ on~$D$ is attained at the
point~$(1/2,t_0)$. In particular, the ground state energy per site is given by
\begin{eqnarray}
  u=u(1/2,t_0)&=\int_0^{x_0(h)}\varepsilon(x)\,\mathrm{d} x-hx_0(h)+\int_0^{1/2}\varepsilon(x)\,\mathrm{d} x-\mu\nonumber\\
  &= \frac1{6}\,\bigg[1-3(h+2\mu)-\frac12(1-4h)^{3/2}\bigg]\,,
      \label{uminB1}
\end{eqnarray}
while the magnetization and charge per site at zero temperature read
\[
  m_s=x_0(h)\,,\qquad n_c=1\,.
\]
The $T=0$ magnetic susceptibility is again
\[
  \chi_s=(1-4h)^{-1/2}=(1-2m_s)^{-1}\,,
\]
while the charge susceptibility vanishes. Since~$s=1/2$ is equivalent to~$N_1+N_2=N$, this is
an~$\mathrm{su}(2)$ phase characterized by the absence of holes.

\medskip
\noindent v)\enspace $0< h/2+\mu<1/4$, $\mu\le h/2$

\smallskip\noindent The above inequalities imply that $u(s,s)=2f_{h/2+\mu}(s)$ attains its unique
global minimum over the interval~$[0,1/2]$ at the point~$s_0=x_0(h/2+\mu)\in(0,1/2)$. It is
straightforward to check that~$u(s,t)$ achieves its minimum over the domain~$D$ at the
point~$(s_0,s_0)$. Indeed, as $h>0$ in this case this minimum cannot be reached on the interior of
the side~$t=0$. Consider next the side~$s=1/2$. If~$h\ge1/4$ the function~$f_h(t)$ is decreasing,
so that
\[
  \fl
  u(1/2,t)=f_h(t)+f_{2\mu}(1/2)\ge f_h(1/2)+f_{2\mu}(1/2)=u(1/2,1/2)>u(s_0,s_0).
\]
On the other hand, if $0\le h<1/4$ then~$f_h$ has a global minimum on~$[0,1/2]$
at~$t_0=x_0(h)\ge x_0(h/2+\mu)=s_0$, with~$t_0\in[0,1/2)$. Since $\varepsilon(x)-2\mu>0$ for~$x> t_0$
(this is obvious for~$\mu\le0$, while for~$0<\mu\le1/8$ it is a consequence of the
inequality~$x_0(2\mu)\le x_0(h)$, which in turn follows from $2\mu\le h$) we have
\[
  f_{2\mu}(1/2)=\int_{0}^{\mathpalette\mathrlapinternal{1/2}}\big(\varepsilon(x)-2\mu\big)\,\mathrm{d} x>
  \int_{0}^{\mathpalette\mathrlapinternal{t_0}}\big(\varepsilon(x)-2\mu\big)\,\mathrm{d} x=f_{2\mu}(t_0)\,.
\]
\begin{table}[t]
  \centering\small
  \begin{tabular}{|c|c|c|}\hline
    Region&
            Ground state energy per site& Spin content\\ \hline\hline
       $\frac h2<\mu<\frac 18$& $\frac1{6}\big[1-3(h+2\mu)-\frac12(1-4h)^{3/2}
                                -\frac12(1-8\mu)^{3/2}\big]$&$B,F_1,F_2$\vru\\ \hline
    $h\ge\frac 14$,\quad $\mu+\frac h2\ge\frac 14$
          & $\frac1{6}\,\big[1-3(h+2\mu)\big]$&$F_1$\vru\\ \hline
        $\mu+\frac h2\le0$& $0$&$B$\vru\\ \hline
        $h< \frac 14$,\quad $\mu\ge\frac 18$&$\frac1{6}\big[1-3(h+2\mu)-\frac12(1-4h)^{3/2}\big]$
                                          &$F_1,F_2$\vru\\ \hline
$\,0<\mu+\frac h2<\frac 14$,\quad $\mu\le\frac h2\,$
          & $\frac16\Big[1-3(h+2\mu)-(1-2h-4\mu)^{3/2}\Big]$&$B,F_1$\vru\\ \hline
  \end{tabular}
  \caption{Ground state energy per site and spin content ($B\equiv$ hole (boson),
    $F_\sigma\equiv$ fermion of type~$\sigma$) for the spin~$1/2$ KY model as a function of its
    parameters~$h\ge0$ and $\mu$ (the unit of energy has been taken as
    $2\pi^2t$).}\label{table.regs}
\end{table}
Hence
\begin{eqnarray*}
  \fl
  u(s_0,s_0)-u(1/2,t_0)&=2f_{h/2+\mu}(s_0)-f_h(t_0)-f_{2\mu}(1/2)\\&<
  2f_{h/2+\mu}(s_0)-f_h(t_0)-f_{2\mu}(t_0)=2f_{h/2+\mu}(s_0)-2f_{h/2+\mu}(t_0)\le0\,,
\end{eqnarray*}
which completes the proof of our assertion. In summary, in this case the ground state energy per
site is given by
\begin{eqnarray}
  u=u(s_0,s_0)&=2\int_0^{\mathpalette\mathrlapinternal{x_0(h/2+\mu)}}
  \varepsilon(x)\,\mathrm{d} x-(h+2\mu)x_0\bigl(\tfrac h2+\mu\bigr)\nonumber\\
  &=\frac16\Big[1-3(h+2\mu)-(1-2h-4\mu)^{3/2}\Big],
  \label{uminB0}
\end{eqnarray}
while the magnetization, the charge density and their susceptibilities (per site) read
\begin{eqnarray*}
  &n_c=2m_s=2x_0(h/2+\mu)=1-\sqrt{1-2h-4\mu}\,,\\
  &\chi_c=4\chi_s=2(1-2h-4\mu)^{-1/2}=\frac2{1-n_c}=\frac{2}{1-2m_s}\,.
\end{eqnarray*}
Note, finally, that the equality~$s=t$ is equivalent to~$N_2=0$. This is thus an~$\mathrm{su}(1|1)$ phase,
consisting only of holes and up fermions.

Our results are summarized in Table~\ref{table.regs}. Note, in particular, that $u$ is continuous
(indeed, of class $C^1$) over its domain, although its second derivatives are discontinuous along
the boundaries of the regions listed in Table~\ref{table.regs} (cf.~Fig.~\ref{fig.u}).

\begin{figure}[t]
  \centering
  \includegraphics[width=.46\textwidth]{hmu.pdf}\hfill
  \includegraphics[width=.46\textwidth]{u.pdf}\hfill
  \caption{Left: phase diagram of the ground state of the spin~$1/2$ KY model ($B\equiv$ hole
    (boson), $F_\sigma\equiv$ fermion of type~$\sigma$). Right: ground state energy per site of the
    latter model as a function of its parameters~$h\ge0$ and $\mu$ (the blue lines represent the
    boundaries of the regions in the phase diagram). In both plots, the unit of energy has been
    taken as~$2\pi^2t$.}
  \label{fig.u}
\end{figure}

\section{Conclusions}\label{sec.conc}

Although the supersymmetric $\mathrm{su}(m)$ $t$-$J$ model with long-range interactions, also known as the
$\mathrm{su}(m)$ KY model, has been studied extensively during the last few decades, an analytical
derivation of its spectrum and exact partition function has been missing so far. With the purpose
of filling up this gap, in this paper we first establish the precise equivalence of the $\mathrm{su}(m)$
KY model to a suitable modification of the $\mathrm{su}(1|m)$ HS spin chain with chemical potential terms.
This equivalence allows us to obtain the partition function of the former model from that of the
latter, which was recently computed in Ref.~\cite{FGLR18}. A remarkable property of this partition
function is the fact that it can be rewritten as the partition function of a suitable
inhomogeneous vertex model.
Analyzing the structure of these two equivalent partition functions, we not only obtain the
complete spectrum of the $\mathrm{su}(m)$ KY model in the presence of an arbitrary magnetic field and
charge chemical potential, but also develop a novel method for determining the magnon numbers or
spin content of the corresponding wave functions. This yields an exhaustive description of the
spectrum in the subspaces with well-defined magnon content in terms of suitably restricted bond
configurations of the equivalent vertex model, which are closely connected with supersymmetric
versions of Haldane's motifs and their related skew Young tableaux. For the particular case $m=2$,
this description provides a rigorous proof of a long-standing conjecture by Saiga and
Kuramoto~\cite{SK99} based on numerical evidence.

As a concrete application of our results, we study various thermodynamic properties of the
$\mathrm{su}(2)$ KY model in the zero temperature limit. To this end, we determine the structure of the
motifs and bond configurations yielding the ground state of the latter model in the thermodynamic
limit for different values of the external parameters. The structure of such bond configurations
leads to a complete description of the different ground state phases in terms of the magnetic
field strength and the charge chemical potential. These phases are characterized by the spin
content (magnon numbers) of the corresponding wave functions, namely an $\mathrm{su}(1|2)$ phase where
holes and fermions of both types co-exist, an $\mathrm{su}(1|1)$ phase with holes and fermions of only one
type, and an $\mathrm{su}(0|2)$ phase with fermions of both types, apart from the trivial phases
consisting of only holes or fermions of one type. We also compute the zero-temperature values of
the energy, magnetization and charge density, along with the magnetic and charge susceptibilities,
for each ground state phase. This description of the thermodynamic properties at zero temperature
goes beyond previously known results, which were derived by different methods and restricted to
the $\mathrm{su}(1|2)$ phase. In particular, our analysis confirms that the strong spin-charge separation
characteristic of the long-range $t$-$J$ model at low temperatures~\cite{KK95} occurs in all
nontrivial phases.

Note, finally, that the description of the spectrum of the~$\mathrm{su}(m)$ KY model in terms of
supersymmetric motifs and their associated Young tableaux derived in this paper makes it possible
to compute in closed form the model's thermodynamic functions at finite temperature, by means of
the transfer matrix method developed in Refs.~\cite{EFG12,FGLR18}. In fact, work on this problem
is currently in progress and shall be presented in a forthcoming publication.

\section*{Acknowledgments}
We would like to thank the anonymous referee, whose remarks helped to improve the presentation.
This work was partially supported by Spain's MINECO grant~FIS2015-63966-P. JAC would also like to
acknowledge the financial support of the Universidad Complutense de Madrid through a 2015
predoctoral scholarship.