# Thermodynamics and criticality of $\mathrm{su}(m)$ spin chains of
Haldane–Shastry type

Federico Finkel<EMAIL_ADDRESS>Artemio González-López<EMAIL_ADDRESS>Depto.
de Física Teórica, Facultad de Ciencias Físicas,
Universidad Complutense de Madrid,
Plaza de las Ciencias 1, 28040 Madrid, SPAIN

###### Abstract

We study the thermodynamics and critical behavior of $\mathrm{su}(m)$ spin
chains of Haldane–Shastry type at zero chemical potential, both in the
$A_{N-1}$ and $BC_{N}$ cases. We evaluate in closed form the free energy per
spin for arbitrary values of $m$, from which we derive explicit formulas for
the energy, entropy and specific heat per spin. In particular, we find that
the specific heat features a single Schottky peak, whose temperature is well
approximated for $m\lesssim 10$ by the corresponding temperature for an
$m$-level system with uniformly spaced levels. We show that at low
temperatures the free energy per spin of the models under study behaves as
that of a one-dimensional conformal field theory with central charge $c=m-1$
(with the only exception of the Frahm–Inozemtsev chain with zero value of its
parameter). However, from a detailed study of the ground state degeneracy and
the low-energy excitations, we conclude that these models are only critical in
the antiferromagnetic case, with a few exceptions that we fully specify.

##  I Introduction

In this paper we shall consider a broad class of spin chains with long-range
interactions modeled on the Haldane–Shastry chain [1, 2], whose interactions
can be expressed in terms of the generators of the $\mathrm{su}(m)$ algebra in
the fundamental representation. More precisely, if $N$ denotes the number of
sites and $m$ is the number of internal degrees of freedom the canonical basis
of the system’s Hilbert space $\otimes_{i=1}^{N}{\mathbb{C}}^{m}$ is spanned
by the vectors

$|s_{1},\dots,s_{N}\rangle:=|s_{1}\rangle\otimes\cdots\otimes|s_{N}\rangle,$
(1)

where $s_{i}\in\\{1,\dots,m\\}$. We define the permutation and spin flip
operators $S_{ij}$ and $S_{i}$ ($1\leqslant i<j\leqslant N$) by the usual
formulas

$\displaystyle S_{ij}|\cdots s_{i}\cdots s_{j}\cdots\rangle$
$\displaystyle=|\cdots s_{j}\cdots s_{i}\cdots\rangle,$ (2) $\displaystyle
S_{i}|\cdots s_{i}\cdots\rangle$ $\displaystyle=|\cdots
m-s_{i}+1\cdots\rangle.$ (3)

The latter operators can be easily expressed in terms of the (Hermitian)
$\mathrm{su}(m)$ generators $T^{\alpha}_{k}$ in the fundamental
representation, where $1\leqslant\alpha\leqslant m^{2}-1$ and the subindex
labels the chain sites. Indeed, using the normalization
$\operatorname{tr}(T_{k}^{\alpha}T_{k}^{\beta})=2\delta_{\alpha\beta}$ we have

$S_{ij}=\frac{1}{m}+\frac{1}{2}\sum_{\alpha=1}^{m^{2}-1}T_{i}^{\alpha}T_{j}^{\alpha}\,,\quad
S_{i}=T_{i}^{\gamma}+\frac{1}{2m}(1-(-1)^{m}),$

where the index $\gamma$ is fixed but arbitrary. In particular, for $m=2$ we
have $T_{k}^{\alpha}=\sigma_{k}^{\alpha}$ and hence 111When $m=2$ we take
$\gamma=1$ in the equation for $S_{i}$ to agree with the standard definition
of the Pauli matrices. $S_{i}=\sigma_{i}^{1}$,
$S_{ij}=\frac{1}{2}(1+{\boldsymbol{\sigma}}_{i}\cdot{\boldsymbol{\sigma}}_{j})$,
where $\sigma^{\alpha}$ ($\alpha=1,2,3$) are the Pauli matrices. Note that the
operators $S_{ij}$ obey the standard permutation algebra

$S_{ij}O_{j}=O_{i}S_{ij},\quad S_{ij}O_{i}=O_{j}S_{ij},\quad
S_{ij}O_{k}=O_{k}S_{ij},$

where $k\neq i,j$ and $O_{k}$ is any operator acting on the $k$-th site.

The first class of spin chains we shall be interested in are the
Haldane–Shastry (HS), Polychronakos–Frahm (PF) [4, 5, 6] and Frahm–Inozemtsev
(FI) [7] chains. They can be collectively defined through the formula 222In
what follows sums over Latin indices will implicitly range form $1$ to $N$,
unless otherwise stated.

$H_{\pm}=\sum_{i<j}J_{ij}(1\mp S_{ij})\,,$ (4)

where

$\displaystyle J_{ij}$
$\displaystyle=\frac{J}{2N^{2}\sin^{2}(\xi_{i}-\xi_{j})},\quad\xi_{k}=\frac{k\pi}{N}\enspace(\text{HS
chain})\,,$ (5) $\displaystyle J_{ij}$
$\displaystyle=\frac{J}{N(\xi_{i}-\xi_{j})^{2}}\,,\quad
H_{N}(\xi_{k})=0\quad(\text{PF chain})\,,$ (6) $\displaystyle J_{ij}$
$\displaystyle=\frac{J}{2N^{2}\sinh^{2}(\xi_{i}-\xi_{j})}\,,\quad
L_{N}^{c-1}(e^{2\xi_{k}})=0\quad(\text{FI chain}).$ (7)

Here $J>0$ is a real constant fixing the energy scale, $H_{N}$ denotes the
Hermite polynomial of degree $N$ and $L_{N}^{c-1}$ is a generalized Laguerre
polynomial of degree $N$ with positive parameter $c$. For reasons that will
become clear in the sequel, we shall sometimes refer to the ferromagnetic
models $H_{+}$ as _bosonic_ and to the antiferromagnetic ones $H_{-}$ as
_fermionic_. From the previous expressions we see that the HS chain can be
naturally considered as a circular chain with equally spaced sites and spin-
spin interactions inversely proportional to the square of the chord distance,
while the PF and FI chains are better regarded as linear chains with sites
$\xi_{k}$ defined in Eqs. (6)-(7) and respectively rational (inverse square)
or hyperbolic interactions.

The spin chains discussed above are all naturally related to the $A_{N-1}$
classical root system. The second type of chain we shall deal with is the
variant of the HS chain related to the $BC_{N}$ root system (HS-B chain),
whose Hamiltonian shall be taken as [9, 10]

$H_{\pm}=\frac{J}{4N^{2}}\sum_{i<j}\bigg{(}\frac{1\mp
S_{ij}}{\sin^{2}(\xi_{i}-\xi_{j})}+\frac{1\mp\widetilde{S}_{ij}}{\sin^{2}(\xi_{i}+\xi_{j})}\bigg{)}\\\
+\frac{J}{8N^{2}}\sum_{i}\left(\frac{\beta_{1}}{\sin^{2}\xi_{i}}+\frac{\beta_{2}}{\cos^{2}\xi_{i}}\right)\big{(}1+\varepsilon
S_{i}\big{)}$ (8a) with $P_{N}^{\beta_{1}-1,\beta_{2}-1}(\cos 2\xi_{k})=0.$
(8b)

Here $J>0$, $P_{N}^{\beta_{1}-1,\beta_{2}-1}$ is a Jacobi polynomial of degree
$N$ and parameters $\beta_{1,2}>0$, $\varepsilon=\pm 1$, and we have used the
abbreviation $\widetilde{S}_{ij}:=S_{ij}S_{i}S_{j}=S_{i}S_{j}S_{ij}$. This
model can be regarded as the _open_ version of the HS chain, with chain sites
$z_{j}:=e^{2\mathrm{i}\xi_{j}}$ lying on the upper unit circle (although in
general not uniformly spaced). The spin at $z_{j}$ interacts both with the
remaining spins at $z_{k}$ (with $k\neq j$) and their reflections
$\bar{z}_{k}$ with respect to the real axis, the interaction strength being
inversely proportional to the square of the (chord) distances $|z_{j}-z_{k}|$
and $|z_{j}-\bar{z}_{k}|$.

The HS spin chain was originally introduced as a parent Hamiltonian for the
Gutzwiller variational ground state for the one-dimensional Hubbard model in
the limit of infinite on-site energy. In fact, the (exact) ground state of the
HS chain contains a Jastrow factor reminiscent of the Laughlin wave function
in the fractional quantum Hall effect [11, 12, 13]. The remarkable properties
of the HS chain became apparent shortly after its introduction. For instance,
it was soon found that its spinon excitations behave as effective particles
obeying Haldane’s fractional statistics [14, 12, 15]. This can indeed be
regarded as the simplest realization of anyons in one dimension. The HS chain
is completely integrable [16, 17, 18], and is actually invariant under the
Yangian quantum algebra $Y(\mathrm{sl}(m))$ even for a finite number of sites,
which in part explains the high degeneracies of its spectrum [19]. It is
closely related to the Wess–Zumino–Novikov–Witten (WZNW) model at level 1
[14], and can also be embedded into a larger class of models constructed from
chiral vertex operators of an appropriate conformal field theory [20, 21].

A characteristic property of the HS chain that distinguishes it from short-
range chains like the Heisenberg model is the fact that it can be obtained as
the strong interaction (large coupling) limit of an integrable one-dimensional
system, namely the (spin) Sutherland model [22, 23, 24]. In fact, the PF and
FI chains can be analogously derived from the integrable spin Calogero [25,
26] and Frahm–Inozemtsev [27] dynamical models. Similarly, the HS-B chain is
the large coupling limit of the spin Sutherland model of $BC_{N}$ type [28, 9,
10]. As first pointed out by Polychronakos [6], the connection between the
spin chains (4)–(8) and the dynamical spin models of Calogero–Sutherland type
mentioned above can be exploited to derive the chains’ partition functions in
closed form [29, 10, 30, 31]. From the common structure of these partition
functions and their relation to the representation theory of the Yangian
algebra in terms of border strips and their associated motifs [32], a
remarkable equivalence between the $A_{N-1}$ chains (4)–(7) and certain
(inhomogeneous) vertex models was established in Ref. [33]. More precisely,
the spectrum of the latter chains (with the correct degeneracy for each
energy) is the same as that of a vertex model with $N+1$ vertices connected by
$N$ bonds, each of which can take the values $1,\dots,m$. The energy of a
configuration of this model, represented by a _bond vector_

$\mathbf{s}:=(s_{1},\dots,s_{N})\in\\{1,2,\dots,m\\}^{N},$

can be computed through the formula

$E(\mathbf{s})=J\sum_{i=1}^{N-1}\delta_{\pm}(s_{i},s_{i+1}){\mathcal{E}}(x_{i})\,,\quad
x_{i}:=i/N\,.$ (9)

Here the _dispersion function_ ${\mathcal{E}}(x)$, which depends on the chain
considered, is given by

${\mathcal{E}}(x)=\begin{cases}x(1-x),\quad&(\text{HS chain})\\\ x,&(\text{PF
chain})\\\ x(x+\gamma_{N}),&(\text{FI chain})\end{cases}$ (10)

with $\gamma_{N}:=(c-1)/N$, while the function $\delta_{\pm}$ (where the $\pm$
sign corresponds to the double sign in $H_{\pm}$) is defined by

$\delta_{+}(i,j)=\begin{cases}0,\enspace s_{i}\leqslant s_{i+1}\\\ 1,\enspace
s_{i}>s_{i+1},\end{cases}\quad\delta_{-}(i,j)=\begin{cases}0,\enspace
s_{i}<s_{i+1}\\\ 1,\enspace s_{i}\geqslant s_{i+1}.\end{cases}$ (11)

We thus see that the spin degrees of freedom behave as bosons (resp. fermions)
in the ferromagnetic (resp. antiferromagnetic) case. Note also that the
vectors with components $\delta_{\pm}(s_{i},s_{i+1})$ are closely connected to
the _motifs_ introduced by Haldane [14]. An analogous description of the
spectrum of the HS-B chain was recently found in Ref. [34]. More precisely, in
this case the vertex model has an additional vertex and a last bond $s_{N+1}$
assuming the _fixed_ (half-integer) value $m_{\varepsilon}+\frac{1}{2}$, where

$m_{\varepsilon}:=\frac{1}{2}(m+\varepsilon\pi(m))$

and $\pi(m)=\frac{1}{2}\,\big{(}1-(-1)^{m}\big{)}$ is the parity of $m$. Thus
in this case

$E(\mathbf{s})=J\sum_{i=1}^{N}\delta_{\pm}(s_{i},s_{i+1}){\mathcal{E}}(x_{i})\,,$
(12)

with $\delta_{\pm}$ as above and dispersion function

${\mathcal{E}}(x)=x\,\bigg{(}\gamma_{N}+1-\frac{x}{2}\bigg{)},\quad\gamma_{N}:=\frac{1}{2N}(\beta_{1}+\beta_{2}-1)\,.$
(13)

The thermodynamics of spin chains of HS type has been studied ever since the
early work of Haldane, who derived an expression for the entropy of the
$\mathrm{su}(2)$ HS chain by means of the spinon description of its spectrum
[35]. Shortly afterwards, Sutherland and Shastry [36] addressed the general
$\mathrm{su}(m)$ case, outlining a complicated procedure for computing the
free energy which involves two successive integrations. This procedure,
however, only yields an explicit expression for $m=2$. Around the same time, a
heuristic formula for the free energy per spin of the PF and FI chains (with
no magnetic field or chemical potential term) was presented in Refs. [5, 7],
again only for the $\mathrm{su}(2)$ case. In fact, a systematic study of the
thermodynamics of the $A_{N-1}$ chains (4)-(7) (with an additional chemical
potential or magnetic field term) using the transfer matrix method was
undertaken for the first time in Ref. [37], and extended later to the
supersymmetric case in Ref. [38]. The key idea in this respect is the fact
that using Eq. (9) (or, actually, its generalization to allow for a chemical
potential term) it is straightforward to express the partition function as the
trace of a product of $N$ site-dependent $m\times m$ transfer matrices. In our
case, it follows from (9) that the partition function of the three $A_{N-1}$
chains can be collectively written as

$Z=\operatorname{tr}\big{[}A(x_{0})A(x_{1})\cdots A(x_{N-1})\big{]}\,,$

where the $m\times m$ _transfer matrix_ $A(x)$ has entries

$A_{\mu\nu}=e^{-\beta J{\mathcal{E}}(x)\delta_{\pm}(\mu,\nu)}\,,\quad
1\leqslant\mu,\nu\leqslant m.$ (14)

Here $\beta:=1/T$ is the inverse temperature (taking Boltzmann’s constant
$k_{\mathrm{B}}$ as unity), and the dispersion relation ${\mathcal{E}}(x)$ is
given by (10). Since the matrix $A(x)$ has positive entries for all
$x\in[0,1]$, the classical Perron theorem [39, 40] implies that $A(x)$ has a
positive _simple_ eigenvalue $\lambda_{1}(x)$ which exceeds the modulus of any
other eigenvalue. From this fact it readily follows that in the thermodynamic
limit $N\to\infty$ the free energy per spin $f(T)$ of the chains (4)-(7) can
be expressed as

$f(T)=-T\int_{0}^{1}\ln\lambda_{1}(x)\,\mathrm{d}x\,.$ (15)

As shown in Ref. [41], similar expressions are valid for the HS-B chain (8).
Indeed, from (12) we obtain

$Z=\operatorname{tr}\big{[}A(x_{1})\cdots A(x_{N-1})B\big{]},$

where $A_{\mu\nu}(x)$ is defined as above but using Eq. (13) for the
dispersion relation, and the $m\times m$ matrix $B$ has entries

$B_{\mu\nu}=e^{-\beta
J(\gamma_{N}+1/2)\delta_{\pm}(\mu,m_{\varepsilon}+\frac{1}{2})}\,,\quad
1\leqslant\mu,\nu\leqslant m\,.$

Since all the matrices in the expression for the partition function have again
positive entries, it follows from Perron’s theorem that the thermodynamic free
energy per spin is given by Eq. (15) also in this case (see [41] for details).
Note that when computing the thermodynamic free energy from Eq. (15) we must
replace the parameter $\gamma_{N}$ in (10)-(13) by

$\gamma:=\lim_{N\to\infty}\gamma_{N}\geqslant 0\,.$

In fact, Eq. (15) has been shown to hold for the $\mathrm{su}(m|n)$
supersymmetric version of the chains (4)–(8) studied in this paper, even with
the addition of a general chemical potential term [38, 41]. Thus the
thermodynamic functions of all of these models can be computed in closed form
provided that the Perron eigenvalue $\lambda_{1}(x)$ of the transfer matrix
$A(x)$ in Eq. (14) can be explicitly found. So far, this has only been done in
the $\mathrm{su}(2)$ case (bosonic or fermionic) [37, 41] and in the
supersymmetric case 333In the truly supersymmetric case the transfer matrix
has always a zero eigenvalue, which is doubly degenerate for $m=n=2$. This
makes it straightforward to diagonalize the latter matrix when $1\leqslant
m,n\leqslant 2$. with $1\leqslant m,n\leqslant 2$ [38, 41].

The main aim of this paper is to derive a remarkably simple expression for
$\lambda_{1}(x)$ for all the HS-type $\mathrm{su}(m)$ chains (4)–(8), valid
for _arbitrary_ values of $m$. Thus the thermodynamic functions of these
models can be evaluated in closed form. We stress that such closed-form
expressions had been obtained so far only for the $\mathrm{su}(2)$ case, even
at zero chemical potential.

Another problem we shall address in this work is the study of the critical
behavior of the chains (4)–(8). As is well known, a strong indication of the
critical character of a model is the low-temperature behavior of its free
energy. The reason for this is that at low temperatures the free energy per
unit length of a $(1+1)$-dimensional conformal field theory (CFT) with central
charge $c$ behaves as [43, 44]

$f(T)\simeq f(0)-\frac{\pi cT^{2}}{6v}$ (16)

where $v$ is the Fermi velocity and we are using natural units $\hbar=1$. It
is thus expected that the free energy of a critical system obey the latter
asymptotic formula at sufficiently low temperatures, with $c$ equal to the
central charge of the effective CFT governing the model’s low energy behavior.
Using the explicit expression (15) for the free energy per spin of the chains
(4)–(8), we shall show that Eq. (16) is satisfied for these models with
central charge $c=m-1$ (with the only exception of the FI chain with
$\gamma=0$). This result agrees with the calculation in Ref. [45] for the
supersymmetric PF chain (using a different method), and is consistent with the
fact that the low energy excitations of the (original) HS chain (5) are
governed by the $\mathrm{su}(m)_{1}$ WZNW model [14, 46, 47]. It should be
emphasized, however, that (16) is only a necessary condition for criticality.
Indeed, a CFT —and thus a truly critical model— must have low-energy
excitations with a linear energy-momentum relation and the degeneracy of its
ground state should be finite. Using Eqs. (9)-(10) and (12)-(13), we shall
prove that both of these conditions hold in our case. In this way we shall
show that the chains (4)–(7) (with $\gamma>0$ for the FI chain) are critical
only in the fermionic (antiferromagnetic) case, which is again in agreement
with the results for the HS chain in Ref. [48]. On the other hand, we shall
see that the HS-B chain is critical not only in the fermionic case, but also
in the bosonic (ferromagnetic) one when $m=2$, or $m=3$ and $\varepsilon=-1$.

We shall close this section with a few words on the paper’s organization. In
Section II we recall the formula for the free energy per spin of
$\mathrm{su}(m)$ chains of HS-type in terms of the Perron eigenvalue of a
suitable transfer matrix derived in Refs. [37] and [41]. We then evaluate this
eigenvalue, thus obtaining a simple closed-form expression for the free energy
of the latter models valid for arbitrary values of $m$. This expression is
used in Section III to derive explicit formulas for the energy, entropy and
specific heat per spin. We also study the approximation of these functions by
those of an $m$-level system with uniformly spaced levels, and establish the
existence of a single Schottky peak in the specific heat for all values of
$m$. In Section IV we analyze the critical behavior of the models under study
by first determining the low-temperature behavior of their free energy, and
then examining in detail the ground state degeneracy and the existence of low-
energy excitations with linear energy-momentum relation. The paper ends with a
technical appendix in which we deduce the full asymptotic expansion of the
free energy per spin at low temperatures used in Section IV.

##  II Free energy

In this section we shall evaluate the Perron (“dominant”) eigenvalue
$\lambda_{1}(x)$ of the transfer matrix $A(x)$ given by (14), which as we have
seen determines the free energy per spin of the HS-type chains (4)–(8) through
the integral (15). To this end, we first recall that according to Perron’s
theorem the dominant eigenvalue of a positive matrix 444A matrix (or, in
particular, a vector) is said to be positive if all its entries entries are
positive. possesses a positive eigenvector. In fact, a corollary of Perron’s
theorem posits that the dominant eigenvalue of a positive matrix is the only
eigenvalue possessing a positive eigenvector. Since this property shall be
essential in what follows, we shall briefly summarize its proof. Indeed,
suppose that $\mathbf{v}$ is a positive eigenvector of a positive matrix $P$
with eigenvalue $\lambda$, and denote by $\lambda_{1}$ the dominant eigenvalue
of $P$. Since the transpose matrix $P^{\mathsf{T}}$ is also positive and has
the same spectrum as $P$, its dominant eigenvalue is also $\lambda_{1}$, and
therefore there exists a positive vector $\mathbf{w}$ such that
$P^{\mathsf{T}}\mathbf{w}=\lambda_{1}\mathbf{w}$. We then have

$\lambda_{1}(\mathbf{w},\mathbf{v})=(P^{\mathsf{T}}\mathbf{w},\mathbf{v})=(\mathbf{w},P\mathbf{v})=\lambda(\mathbf{w},\mathbf{v}),$

where we have used the fact that $P^{\mathsf{T}}$ and $\lambda_{1}$ are real.
Since both $\mathbf{v}$ and $\mathbf{w}$ are positive vectors, from the latter
equality it follows that $\lambda=\lambda_{1}$, as claimed.

In fact, if an eigenvalue $\lambda$ of a positive matrix $P$ is known,
Perron’s theorem provides a simple test to ascertain whether $\lambda$ is the
dominant eigenvalue of $P$. Indeed, if we denote by $C^{\lambda}_{ij}$ the
$(i,j)$ cofactor of the matrix $\lambda-P$ we then have the elementary
identity

$\sum_{j=1}^{m}(\lambda\delta_{ij}-P_{ij})C^{\lambda}_{kj}=\delta_{ik}\det(\lambda-P)=0,$

where $m$ is the order of $P$. From the previous identity it follows that
_any_ row of the cofactor matrix $(C_{ij}^{\lambda})_{i,j=1}^{m}$ is an
eigenvector of $P$ with eigenvalue $\lambda$. However, as part of the proof of
Perron’s theorem, it is shown that _all_ the cofactors $C_{ij}^{\lambda_{1}}$
corresponding to the dominant eigenvalue $\lambda_{1}$ are _positive_ [40].
From the discussion in the previous paragraph we then obtain the following
elementary test: an eigenvalue $\lambda$ of a positive matrix $P$ is its
dominant eigenvalue if and only if _any_ row of the cofactor matrix of
$\lambda-P$ is a positive vector.

Let us now turn to the computation of the dominant eigenvalue of the positive
matrix $A(x)$ in Eq. (14). To begin with, from the definition (11) of
$\delta_{+}$ it follows that in the _bosonic_ case the matrix $A(x)$ has the
following structure

$A(x)=\begin{pmatrix}1&1&\cdots&1&1\\\ a^{m}&1&\cdots&1&1\\\
\vdots&\vdots&\ddots&\vdots&\vdots\\\ a^{m}&a^{m}&\cdots&1&1\\\
a^{m}&a^{m}&\cdots&a^{m}&1\end{pmatrix},$

where we have set 555For the sake of conciseness, we have omitted the
dependence of $a$ on $x$.

$a:=e^{-\beta J{\mathcal{E}}(x)/m}\,.$

An elementary calculation shows that

$\lambda(x)=\sum_{k=0}^{m-1}a^{k}=\frac{1-a^{m}}{1-a}$

is an eigenvalue of $A(x)$. We could now apply the previous test to check
whether $\lambda(x)$ is the dominant eigenvalue of $A(x)$, but in this case it
is easier to observe that the vector

$\mathbf{v}=(1,a,\dots,a^{m-1})$

is a positive eigenvector of $A(x)$ with eigenvalue $\lambda(x)$. Indeed,

$\displaystyle[A(x)\mathbf{v}]_{k}$
$\displaystyle=\sum_{i=1}^{k-1}a^{m}a^{i-1}+\sum_{i=k}^{m}a^{i-1}=\sum_{i=k-1}^{m+k-2}a^{i}$
$\displaystyle=a^{k-1}\lambda(x)=\lambda(x)v_{k}\,.$

From the previous discussion it follows that $\lambda(x)=\lambda_{1}(x)$. We
thus obtain the following remarkable formula for the free energy per spin
$f_{+}(T)$ of the _bosonic_ HS-type chains (4)–(8):

$\displaystyle f_{+}(T)$
$\displaystyle=-T\int_{0}^{1}\ln\left(\sum_{k=0}^{m-1}e^{-\beta
J{\mathcal{E}}(x)k/m}\right)\\!\mathrm{d}x$
$\displaystyle=-T\int_{0}^{1}\ln\biggl{(}\frac{1-e^{-\beta
J{\mathcal{E}}(x)}}{1-e^{-\beta J{\mathcal{E}}(x)/m}}\biggr{)}\mathrm{d}x$
$\displaystyle=f_{0}-T\int_{0}^{1}\ln\left[\frac{\sinh\Bigl{(}\beta
J{\mathcal{E}}(x)/2\Bigr{)}}{\sinh\Bigl{(}\beta
J{\mathcal{E}}(x)/(2m)\Bigr{)}}\right]\\!\mathrm{d}x,$ (17)

where

$f_{0}:=\frac{J}{2}\bigg{(}1-\frac{1}{m}\bigg{)}{\mathcal{E}}_{0}\,,\quad{\mathcal{E}}_{0}:=\int_{0}^{1}{\mathcal{E}}(x)\mathrm{d}x\,.$

The constant ${\mathcal{E}}_{0}$, which can be easily found from (10)-(13),
takes the values $1/6$ (HS chain), $1/2$ (PF chain), and $\gamma/2+1/3$ (FI
and HS-B chains).

An analogous result for the fermionic case can be easily derived noting that

$A_{-}(x)=A_{+}(x)+(a^{m}-1)\mathbb{1}\,,$

where for clarity’s sake we have denoted by $A_{-}(x)$ (resp. $A_{+}(x)$) the
matrix $A(x)$ in the fermionic (resp. bosonic) case. Hence the dominant
eigenvalue of the matrix $A_{-}(x)$ is given by

$\lambda_{1}(x)=\sum_{k=0}^{m-1}a^{k}+a^{m}-1=a\sum_{k=0}^{m-1}a^{k}\,.$

Substituting into (15) we obtain the following simple formula for the free
energy per spin $f_{-}(T)$ in the _fermionic_ case:

$f_{-}(T)=f_{+}(T)+\frac{J{\mathcal{E}}_{0}}{m}\,.$ (18)

For $m=2$, the previous formulas for $f_{\pm}$ coincide with those in Refs.
[37, 41]. As an additional consistency check, note that as $T\to\infty$ from
Eqs. (17)-(18) we easily obtain

$f(T)\simeq-T\ln m\,,$

in agreement with the elementary identity

$Z(T)\underset{T\to\infty}{\simeq}m^{N}\,.$

######  Remark 1.

The thermodynamic free energy of the PF chain can be computed in closed form
in terms of the dilogarithm function [51]

$\operatorname{Li}_{2}(z):=-\int_{0}^{z}\frac{\ln(1-t)}{t}\,\mathrm{d}t\,,$

analytic in the cut complex plane ${\mathbb{C}}\setminus[1,\infty)$, where the
integral is taken along any path in the latter set joining the origin to the
point $z$. Indeed, we have

$f_{+}(T)=\frac{T^{2}}{J}\bigg{[}m\operatorname{Li}_{2}(e^{-\beta
J/m})-\operatorname{Li}_{2}(e^{-\beta J})-\frac{\pi^{2}}{6}\,(m-1)\bigg{]}.$

##  III Thermodynamics

The energy density $u$, entropy per spin $s$ and specific heat per spin
$c_{\mathrm{V}}$ of the HS-type chains (4)–(8) can be readily computed using
Eqs. (17)-(18). Before doing so, to simplify our formulas we shall set without
loss of generality $J=1$, so that energy and temperature become dimensionless.
With this proviso, we find

$\displaystyle u_{\pm}$ $\displaystyle=\frac{\partial}{\partial\beta}(\beta
f_{\pm})=\frac{1}{2}\bigg{(}1\mp\frac{1}{m}\bigg{)}{\mathcal{E}}_{0}-\frac{1}{2}\int_{0}^{1}{\mathcal{E}}(x)\left[\coth\biggl{(}\frac{\beta{\mathcal{E}}(x)}{2}\biggr{)}-\frac{1}{m}\coth\biggl{(}\frac{\beta{\mathcal{E}}(x)}{2m}\biggr{)}\right]\\!\mathrm{d}x,$
(19) $\displaystyle s$
$\displaystyle=\beta(u_{\pm}-f_{\pm})=\int_{0}^{1}\left\\{\ln\left[\frac{\sinh\left(\frac{\beta{\mathcal{E}}(x)}{2}\right)}{\sinh\left(\frac{\beta{\mathcal{E}}(x)}{2m}\right)}\right]-\frac{\beta{\mathcal{E}}(x)}{2}\bigg{[}\coth\biggl{(}\frac{\beta{\mathcal{E}}(x)}{2}\biggr{)}-\frac{1}{m}\coth\biggl{(}\frac{\beta{\mathcal{E}}(x)}{2m}\biggr{)}\bigg{]}\right\\}\\!\mathrm{d}x,$
(20) $\displaystyle c_{\mathrm{V}}$ $\displaystyle=-\beta^{2}\frac{\partial
u_{\pm}}{\partial\beta}=\frac{\beta^{2}}{4}\int_{0}^{1}{\mathcal{E}}(x)^{2}\left[\frac{1}{m^{2}}\operatorname{csch}^{2}\biggl{(}\frac{\beta{\mathcal{E}}(x)}{2m}\biggr{)}-\operatorname{csch}^{2}\biggl{(}\frac{\beta{\mathcal{E}}(x)}{2}\biggr{)}\right]\\!\mathrm{d}x.$
(21)

Using the explicit formula for the free energy per spin of the PF chain from
the previous section we readily find the following simple expressions for its
thermodynamic functions:

$\displaystyle u_{+}$
$\displaystyle=u_{-}-\frac{1}{2m}=-f_{+}-T\ln\left(\frac{1-e^{-\beta}}{1-e^{-\beta/m}}\right),$
$\displaystyle s$ $\displaystyle=-2\beta
f_{+}-\ln\left(\frac{1-e^{-\beta}}{1-e^{-\beta/m}}\right),$ $\displaystyle
c_{\mathrm{V}}$ $\displaystyle=-2\beta
f_{+}-2\ln\left(\frac{1-e^{-\beta}}{1-e^{-\beta/m}}\right)+\beta\left(\frac{1}{e^{\beta}-1}-\frac{1/m}{e^{\beta/m}-1}\right).$

From the first equality in Eq. (17) it follows that for large $T$ (i.e., for
$T\gg{\mathcal{E}}(1)$ for the PF, FI and HS-B chains or
$T\gg{\mathcal{E}}(1/2)=1/4$ for the HS chain) the ferromagnetic free energy
per spin $f_{+}(T)$ can be well approximated replacing ${\mathcal{E}}(x)$ by
its mean value ${\mathcal{E}}_{0}$ over the interval $[0,1]$. In other words,
we have

$f_{+}(T)\underset{T\to\infty}{\simeq}-T\ln\left(\sum_{k=0}^{m-1}e^{-k\beta\frac{{\mathcal{E}}_{0}}{m}}\right)=:f_{m}(T).$

The right-hand side is the partition function of an $m$-level system with
uniformly spaced levels $E_{k}=k{\mathcal{E}}_{0}/m$, $k=0,\dots,m-1$. In
fact, at sufficiently high temperatures the thermodynamic functions of the HS-
type chains (4)-(8) behave qualitatively as those of the corresponding $m$
level system, as can be seen, for instance, from Fig. 1 for the HS chain with
$m=2,\dots,5$ (we omit the corresponding plots for the PF, FI and HS-B chains,
which are very similar). It should be noted, however, that at low temperatures
the thermodynamic functions of the HS-type chains behave quite differently
than those of their corresponding $m$-level system. Indeed, at low
temperatures the free energy of the $m$-level system,

$\displaystyle f_{m}(T)$
$\displaystyle=-T\ln\left(\frac{1-e^{-\beta{\mathcal{E}}_{0}}}{1-e^{-\beta{\mathcal{E}}_{0}/m}}\right)$
$\displaystyle=-Te^{-\beta{\mathcal{E}}_{0}/m}+O(Te^{-2\beta{\mathcal{E}}_{0}/m}),$

is exponentially small, and so are its remaining thermodynamic functions. On
the other hand, from the discussion in the next section (cf. Eq. (26)) it
follows that for the HS-type chains $f_{\pm}(T)-f_{\pm}(0)\sim-T^{2}$ as $T\to
0+$ (or $\sim-T^{3/2}$ for the FI chain with $\gamma=0$; cf. Eq. (34)). Thus
$u-u(0)\sim T^{2}$, $s\sim T$ and $c_{V}\sim T$ as $T\to 0+$ (or $u-u(0)\sim
T^{3/2}$, $s\sim T^{1/2}$ and $c_{V}\sim T^{1/2}$ for the FI chain with
$\gamma=0$).

Figure 1: Thermodynamic free energy, energy, entropy, and specific heat per
spin of the $\mathrm{su}(m)$ HS chain (4)-(5) with $m=2,3,4,5$ (green, red,
blue, and orange solid lines, respectively), compared to their counterparts
for an $m$-level system (dashed lines).

From Fig. 1 (and similar plots for the PF, FI and HS-B chains that we are not
displaying) it is apparent that, at least for low values of $m$, the specific
heat per spin of the HS-type chains features a single absolute maximum (i.e.,
a Schottky peak). This behavior can be qualitatively understood by analyzing
the specific heat of the corresponding $m$-level system, given by

$c_{\mathrm{V},m}=t^{2}\left[\frac{1}{m^{2}}\operatorname{csch}^{2}(t/m)-\operatorname{csch}^{2}t\right],\quad
t:=\frac{\beta{\mathcal{E}}_{0}}{2}.$

Indeed, although we have not found an analytic proof, we have compelling
numerical evidence that $c_{\mathrm{V},m}(t)$ has a single critical point
$t_{m}$ in the half-line $(0,\infty)$ which is an absolute maximum. It is
therefore to be expected that the specific heat of the $\mathrm{su}(m)$ HS-
type chains feature a single Schottky peak at a temperature $T_{m}$ of the
order of ${\mathcal{E}}_{0}/(2t_{m})$. In fact, for low values of $m\lesssim
10$ the Schottky temperature $T_{m}$ is reasonably well approximated by the
$m$-level value ${\mathcal{E}}_{0}/(2t_{m})$, as can be seen from Fig. 2
(left) for the HS chain (and similarly for the other HS-type chains). On the
other hand, for large $m$ the difference ${\mathcal{E}}_{0}/(2t_{m})-T_{m}$
appears to tend to a constant positive value, while $T_{m}\to 0$ (see Fig. 2,
right).

Figure 2: Left: Temperature $T_{m}$ of the Schottky peak of the
$\mathrm{su}(m)$ HS chain with $m=2,\dots,200$ compared to its rough $m$-level
approximation ${\mathcal{E}}_{0}/(2t_{m})=1/(12t_{m})$. Right: difference
${\mathcal{E}}_{0}/(2t_{m})-T_{m}$ for the HS, PF, FI HS-B chains (with
$\gamma=0$ for the last two chains) for $m=2,\dots,200$.

##  IV Critical behavior

As mentioned in the Introduction, one of the hallmarks of criticality is the
low-temperature behavior of the free energy, given by Eq. (16). In our case,
using Eqs. (17)-(18) with $J=1$ we obtain

$\displaystyle f_{\pm}(T)-f_{\pm}(0)=$
$\displaystyle-T\int_{0}^{1}\ln\left(1-e^{-\beta{\mathcal{E}}(x)}\right)\mathrm{d}x$
$\displaystyle+T\int_{0}^{1}\ln\left(1-e^{-\beta{\mathcal{E}}(x)/m}\right)\mathrm{d}x\,.$
(22)

As $T\to 0+$, the main contribution to both of these integrals comes from a
small neighborhood of the points where the dispersion relation
${\mathcal{E}}(x)$ vanishes, i.e., $x=0,1$ for the HS chain and $x=0$ for the
PF, FI and HS-B chains. For the latter three chains ${\mathcal{E}}(x)$ is
monotonically increasing over the interval $[0,1]$, so that performing the
changes of variables $y=\beta{\mathcal{E}}(x)$ and $y=\beta{\mathcal{E}}(x)/m$
in the integrals in (22) we have

$f_{\pm}(T)-f_{\pm}(0)=-T^{2}\int_{0}^{\beta{\mathcal{E}}(1)}\ln\left(1-e^{-y}\right)({\mathcal{E}}^{-1})^{\prime}(Ty)\,\mathrm{d}y\\\
+mT^{2}\int_{0}^{\frac{\beta{\mathcal{E}}(1)}{m}}\ln\left(1-e^{-y}\right)({\mathcal{E}}^{-1})^{\prime}(mTy)\,\mathrm{d}y,$
(23)

where ${\mathcal{E}}^{-1}:0\to{\mathcal{E}}(1)$ denotes the inverse function
of ${\mathcal{E}}:0\to 1$. Since the main contribution to both integrals comes
from the point $y=0$, we can approximate $({\mathcal{E}}^{-1})^{\prime}(Ty)$
and $({\mathcal{E}}^{-1})^{\prime}(mTy)$ by
$({\mathcal{E}}^{-1})^{\prime}(0)=1/{\mathcal{E}}^{\prime}(0)$ (assuming that
$\gamma>0$ for the FI chain) and push the upper limit in each integral to
$+\infty$, thus obtaining

$\displaystyle f_{\pm}(T)-f_{\pm}(0)$
$\displaystyle=\frac{(m-1)T^{2}}{{\mathcal{E}}^{\prime}(0)}\int_{0}^{\infty}\ln\left(1-e^{-y}\right)\,\mathrm{d}y+O(T^{3})$
$\displaystyle=-\frac{(m-1)T^{2}}{{\mathcal{E}}^{\prime}(0)}\,\operatorname{Li}_{2}(1)+O(T^{3})$
$\displaystyle=-\frac{(m-1)\pi^{2}}{6{\mathcal{E}}^{\prime}(0)}\,T^{2}+O(T^{3});$
(24)

see the appendix for more details 666In fact, it is shown in the appendix that
for the PF chain the $O(T^{3})$ term is actually $O(T^{2}e^{-\beta/m})$.. It
follows from (24) that the free energy per spin of the PF, FI (with
$\gamma>0$) and HS-B chains behaves as that of a CFT with central charge
$c=m-1$. To see this, note first of all that the variable $x$ can be regarded
as $p/\pi$, where $p$ is the momentum (defined modulo $2\pi$). Indeed, since
the dispersion relation is monotonic the interval $0\leqslant x\leqslant 1$
corresponds to the positive momentum range $0\leqslant p\leqslant\pi$ and not
to the full range $-\pi\leqslant p\leqslant\pi$. As the relation between
energy ${\mathcal{E}}(x)$ and momentum $p=\pi x$ is linear near $p=0$, the
Fermi velocity is given by

$v=\left.\frac{\mathrm{d}{\mathcal{E}}}{\mathrm{d}p}\right|_{p=0}=\frac{{\mathcal{E}}^{\prime}(0)}{\pi}\,,$
(25)

so that (24) can indeed be written as

$f_{\pm}(T)-f_{\pm}(0)=-\frac{(m-1)\pi}{6v}\,T^{2}+O(T^{3}).$ (26)

For the HS chain (4)-(5), the dispersion relation ${\mathcal{E}}(x)$ is
symmetric about $x=1/2$ and increasing in the interval $[0,1/2]$. Hence we can
write

$f_{\pm}(T)-f_{\pm}(0)=-2T\int_{0}^{1/2}\ln\left(1-e^{-\beta{\mathcal{E}}(x)}\right)\mathrm{d}x\\\
+2T\int_{0}^{1/2}\ln\left(1-e^{-\beta{\mathcal{E}}(x)/m}\right)\mathrm{d}x\,,$
(27)

and proceeding as above we arrive at

$f_{\pm}(T)-f_{\pm}(0)=-\frac{(m-1)\pi^{2}}{3{\mathcal{E}}^{\prime}(0)}\,T^{2}+O(T^{3})$
(28)

(see again the appendix for details on the error term). However, in this case
the symmetry of the dispersion relation about $x=1/2$ implies that the
relation between $x$ and $p$ is $p=2\pi x$ (the interval $0\leqslant
x\leqslant 1$ now corresponds to the full momentum range $0\leqslant
p\leqslant 2\pi$). Hence

$v=\left.\frac{\mathrm{d}{\mathcal{E}}}{\mathrm{d}p}\right|_{p=0}=\frac{{\mathcal{E}}^{\prime}(0)}{2\pi}\,,$
(29)

and (26) also holds in this case. Thus Eq. (26) is valid for all the HS-type
chains, except for the FI chain with $\gamma=0$. (In fact, as shown in the
appendix, for the latter chain $f_{\pm}(T)-f_{\pm}(0)\sim-T^{3/2}$ as $T\to
0+$.)

As mentioned in the Introduction, to ascertain the criticality of a quantum
system we must also examine the degeneracy of its ground state and study its
low-energy excitations. Both of these problems can be addressed with the help
of Eqs. (9)-(10) and (12)-(13) for the energy spectrum. We shall start by
determining the ground state of each of the chains (4)-(8) and its degeneracy.

Consider, first, the PF, FI and HS-B chains, whose dispersion relation
${\mathcal{E}}(x)$ is monotonically increasing over the interval $[0,1]$. In
the fermionic case $\delta(s_{i},s_{i}+1)=0$ if and only if $s_{i}<s_{i+1}$,
and hence the bond vectors $\mathbf{s}_{g}$ yielding the ground state of the
PF and FI chains are obtained by placing 777We use the notation $\lfloor
x\rfloor$ to denote the integer part of the real number $x$. $r:=\left\lfloor
N/m\right\rfloor$ sequences $(1,\dots,m)$ starting from the right end and
filling the cremaining $N-rm$ components with an increasing sequence
$(s_{1},\dots,s_{N-rm})\in\\{1,\dots,m\\}^{N-rm}$:

$\displaystyle\mathbf{s}_{g}=(s_{1},\dots,s_{N-rm},1,\dots,m,\dots,1,\dots,m),$
(30) $\displaystyle 1<s_{1}<\cdots<s_{N-rm}\,.$

The ground state degeneracy is therefore $\binom{m}{N-rm}\ll N$ as
$N\to\infty$. In the case of the HS-B chain, since
$\delta(s_{N},s_{N+1})=\delta(s_{N},m_{\varepsilon}+\frac{1}{2})$, to obtain
the ground state we must take $s_{N}=m_{\varepsilon}$. Hence the ground state
bond vectors are in this case

$\displaystyle\mathbf{s}_{g}=(s_{1},\dots,s_{N-rm},1,\dots,m,\dots,1,\dots,m,1,\dots,m_{\varepsilon}),$
$\displaystyle 1<s_{1}<\cdots<s_{N-rm},$

where now $r=\left\lfloor(N-m_{\varepsilon})/m\right\rfloor$. The ground state
degeneracy is again $\binom{m}{N-rm}\ll N$. Thus the degeneracy of the
fermionic PF, FI and HS-B chains remains finite in the thermodynamic limit.
The situation is completely different in the bosonic case, since now
$\delta(s,s)=0$. Thus the ground state bond vectors of the PF and FI chains
are of the form

$\mathbf{s}_{g}=(\underbrace{s_{1},\dots,s_{1}}_{k_{1}},\dots,\underbrace{s_{r},\dots,s_{r}}_{k_{r}})$
(31)

with $k_{1}+\cdots+k_{r}=N$ and $1\leqslant s_{1}<\cdots<s_{r}\leqslant m$.
Thus in this case the ground state degeneracy is given by

$d_{g}=\sum_{r=1}^{m}{\mathcal{P}}(N;r)\binom{m}{r},$

where ${\mathcal{P}}(N;r)$ denotes the number of partitions of the integer $N$
in $r$ parts. In particular, the ground state degeneracy clearly tends to
infinity in the thermodynamic limit. On the other hand, for the bosonic HS-B
chain we must impose the additional restriction $s_{r}\leqslant
m_{\varepsilon}$, so that in this case we have

$d_{g}=\sum_{r=1}^{m_{\varepsilon}}{\mathcal{P}}(N;r)\binom{m_{\varepsilon}}{r}.$

It follows that in this case the ground state is non-degenerate if and only if
$m_{\varepsilon}=1$, i.e., for $m=2$ or $m=3$ and $\varepsilon=-1$. The ground
state degeneracy is otherwise infinite (at least $N+1$) in the thermodynamic
limit.

The above analysis must be slightly modified in the case of the HS chain,
whose dispersion relation is increasing over the interval $[0,1/2]$ and
symmetric about $1/2$. For this reason, to the bond vectors (30) we should add
their “reflected” counterparts

$(1,\dots,m,\dots,1,\dots,m,s_{1},\dots,s_{N-rm}).$

Thus in this case the ground state degeneracy is again finite in the
thermodynamic limit ($2\binom{m}{N-rm}$ when $N$ is not a multiple of $m$, or
$1$ when it is). Finally, in the bosonic case the ground state bond vector is
still of the form (31), and hence the ground state degeneracy is infinite in
the thermodynamic limit. In summary, from the analysis of the ground state
degeneracy we conclude that only the fermionic chains (4)-(8), and the bosonic
HS-B chain with $m=2$ or $m=3,\varepsilon=-1$, can be truly critical. It is
interesting to observe in this respect that the bosonic $\mathrm{su}(3)$ HS-B
chain with $\varepsilon=1$ has exactly the same thermodynamic functions as its
counterpart with $\varepsilon=-1$, but only the latter can be critical.

Let us now turn to the study of the low energy excitations over the ground
state in the possible critical cases identified in the previous paragraph. To
begin with, consider the fermionic $\mathrm{su}(m)$ HS chain. A low energy
excitation over a ground state with bond vector $\mathbf{s}_{g}$ is obtained,
for instance, replacing a component $s_{i}$ with $i\ll N$ (or $N-i\ll N$) and
$s_{i}<s_{i+1}$, by $s^{\prime}_{i}\geqslant s_{i+1}$. In the first case we
add an energy

$\Delta E={\mathcal{E}}(x_{i})\simeq{\mathcal{E}}^{\prime}(0)x_{i}=O(1/N),$

while in the second one

$\Delta
E={\mathcal{E}}(x_{i})={\mathcal{E}}(1-x_{i})\simeq{\mathcal{E}}^{\prime}(0)(1-x_{i})=O(1/N).$

On the other hand, it is well known [14] that the state described by a bond
vector $\mathbf{s}$ has momentum

$P(\mathbf{s})=2\pi\sum_{i=1}^{N-1}x_{i}\delta(s_{i},s_{i+1})\quad\mod{2\pi}.$

Hence by replacing $s_{i}$ by $s_{i}^{\prime}$ we add a momentum

$\Delta P=\pm 2\pi x_{i}\quad\mod{2\pi}\,$

where the “$+$” sign corresponds to the case $i\ll N$. Thus in this case the
model is critical, and its Fermi velocity is given by

$v=\lim_{N\to\infty}\left|\frac{\Delta E}{\Delta
P}\right|=\frac{{\mathcal{E}}^{\prime}(0)}{2\pi},$

in agreement with Eq. (29). Note that this conclusion is consistent with the
discussion in Ref. [48] for the $\mathrm{su}(m|n)$-supersymmetric HS chains.

For the fermionic PF, FI and HS-B chains the situation is quite different, as
the Hamiltonian of these models is not invariant under translations along the
lattice and thus momentum is not conserved. Of course, these models still
possess low energy excitations obtained by exciting a component of the ground
state bond vector with index $i\ll N$ (both in the fermionic case and for the
bosonic HS-B chain with $m=2$ or $m=3,\varepsilon=-1$), with energy

$\Delta E={\mathcal{E}}(x_{i})\simeq{\mathcal{E}}^{\prime}(0)x_{i}=O(1/N).$

By analogy with the HS chain, we assign to an energy eigenstate with bond
vector $\mathbf{s}$ an effective momentum

$P(\mathbf{s})=\pi\sum_{i=1}^{N-\eta}x_{i}\delta(s_{i},s_{i+1})\quad\mod{2\pi},$
(32)

where $\eta=0$ for the HS-B chain and $\eta=1$ otherwise. Note that, as
explained above, in this case the factor multiplying the sum is $\pi$ instead
of $2\pi$, since the dispersion relation of the PF, FI, and HS-B chains is
monotonically increasing over the interval $[0,1]$, and thus this interval
represents only the positive momentum range. With this definition the change
in effective momentum of the low energy excitations described above is $\Delta
P=\pi x_{i}$, and hence their Fermi velocity is given by Eq. (25). In summary,
the HS-type chains chains (4)-(8) are all critical in the fermionic case,
while in the bosonic case only the HS-B chain is critical when $m=2$ or
$m=3,\varepsilon=-1$.

######  Remark 2.

The PF chain is known to possess Yangian invariance [54], and it is
conjectured that the same is true for the FI chain in view of the structure of
its partition function and the high degeneracy of its spectrum [31]. Likewise,
the HS-B chain is also known to have twisted Yangian symmetry [9]. It is thus
conceivable that the effective momentum (32) could be related to the
eigenvalues of one of the conserved Yangian generators for these models.

##  V Conclusions

In this paper we have completely determined in closed form the thermodynamics
of the $\mathrm{su}(m)$ spin chains of Haldane–Shastry type (with zero
chemical potential) (4)-(8) for all $m$. Our method relies on the fact that
the energy spectrum of these models coincides with that of an appropriate
vertex model, which makes it possible to express their thermodynamic free
energy as an integral involving the Perron eigenvalue of a position-dependent
$m\times m$ transfer matrix. We have been able to compute this eigenvalue in
closed form for arbitrary values of $m$ by applying the classical Perron
theorem on positive matrices and some of its consequences. This yields an
explicit expression for the free energy per spin of these models, which fully
determines their thermodynamics. We have found that at sufficiently high
temperatures the thermodynamic functions of the $\mathrm{su}(m)$ HS-type
chains behave qualitatively as those of an $m$-level system with uniformly
spaced levels. In particular, for all values of $m$ the specific heat features
a single Schottky peak, whose temperature is close to that of the
corresponding $m$-level system for $m\lesssim 10$.

Using our explicit formula for the free energy per spin, we have also examined
the critical behavior of the $\mathrm{su}(m)$ HS-type chains (4)-(8). We have
first shown that the low-temperature behavior of the free energy coincides
with that of a CFT with central charge $c=m-1$, both in the ferromagnetic
(bosonic) and the antiferromagnetic (fermionic) regimes (with $\gamma>0$ for
the FI chain). This low temperature behavior of the free energy is, however,
only a necessary condition for criticality. To ascertain whether the models
under study are critical or not, we have used the motif-based description of
their spectrum to study the degeneracy of their ground state and the existence
of low-energy excitations with a linear energy-momentum relation. In this way
we have shown that the antiferromagnetic chains are all critical (again, with
$\gamma>0$ for the FI chain), whereas in the ferromagnetic case only the
$\mathrm{su}(2)$ and $\mathrm{su}(3)$ HS-B chain (with $\varepsilon=-1$ in the
latter case) are critical.

Our main result does not appear to be easily generalizable to the case of non-
zero chemical potential, and in fact the known expressions of the
thermodynamic functions in the latter case for low values of $m$ have a more
complicated structure than those found in this paper. Note, however, that our
closed-form expression for the Perron eigenvalue in the zero chemical
potential case could be used to identify and approximate this eigenvalue for
sufficiently small values of the chemical potentials, which would in turn
yield a corresponding approximation for the thermodynamic functions.

## Acknowledgments

This work was partially supported by grants PGC2018-094898-B-I00 from Spain’s
Ministerio de Ciencia, Innovación y Universidades and G/6400100/3000 from
Universidad Complutense de Madrid.

##  Appendix A Asymptotic expansion of the free energy

In this appendix we shall provide the details of the computation of the
asymptotic approximation to the free energy per spin of the chains (4)–(8) at
low temperatures (see, e.g., Eq. (24)). In fact, we shall derive a complete
asymptotic expansion of the latter function for $T\to 0+$.

Consider, to begin with, the FI chain with $\gamma=0$, for which
${\mathcal{E}}(x)=x^{2}$ and hence ${\mathcal{E}}^{-1}(y)=\sqrt{y}$,
$0\leqslant y\leqslant 1$. From Eq. (23) we have

$f_{\pm}(T)-f_{\pm}(0)=-T^{2}\big{(}g(T)-mg(mT)\big{)}$ (33)

with

$g(T)=\frac{1}{2}\,T^{-1/2}\int_{0}^{\beta}\ln(1-e^{-y})\,\frac{\mathrm{d}y}{\sqrt{y}}\,.$

Since

$\left|\int_{0}^{\infty}\ln(1-e^{-y})\,\frac{\mathrm{d}y}{\sqrt{y}}-\int_{0}^{\beta}\ln(1-e^{-y})\,\frac{\mathrm{d}y}{\sqrt{y}}\right|\\\
=-\int_{\beta}^{\infty}\ln(1-e^{-y})\,\frac{\mathrm{d}y}{\sqrt{y}}=O(T^{1/2}e^{-\beta})$

we can write

$\displaystyle g(T)$
$\displaystyle=\frac{1}{2}\,T^{-1/2}\int_{0}^{\infty}\ln(1-e^{-y})\,\frac{\mathrm{d}y}{\sqrt{y}}+O(e^{-\beta})$
$\displaystyle=-\frac{\sqrt{\pi}}{2}\,\zeta(3/2)T^{-1/2}+O(e^{-\beta}),$

where $\zeta(z):=\sum_{n\geqslant 1}n^{-z}$ is Riemann’s zeta function. From
Eq. (33) we finally obtain

$f_{\pm}(T)-f_{\pm}(0)\\\
=-\left(\sqrt{m}-1\right)\frac{\sqrt{\pi}}{2}\,\zeta(3/2)T^{3/2}+O(T^{2}e^{-\beta/m}),$
(34)

which shows that the FI chain is not critical for $\gamma=0$.

Consider next the PF chain. In this case we simply have

$\int_{0}^{1}\ln(1-e^{-\beta
x})\,\mathrm{d}x=T\int_{0}^{\beta}(1-e^{-y})\,\mathrm{d}y=-\frac{\pi^{2}T}{6}+O(Te^{-\beta}),$

and from (22) we find

$f_{\pm}(T)-f_{\pm}(0)=-(m-1)\frac{\pi^{2}T^{2}}{6}+O(T^{2}e^{-\beta/m}).$

This is in agreement with Eq. (24), since in this case ${\mathcal{E}}(x)=x$.

Let us next deal with the FI and HS-B chains with $\gamma>0$, for which we
respectively have

$({\mathcal{E}}^{-1})^{\prime}(z)=\begin{cases}(\gamma^{2}+4z)^{-1/2},&0\leqslant
z\leqslant\gamma+1\\\\[3.0pt]
\left[(\gamma+1)^{2}-2z)\right]^{-1/2},&0\leqslant
z\leqslant\gamma+\frac{1}{2}.\end{cases}$

Since in both cases all the derivatives of ${\mathcal{E}}(z)$ are bounded in
their respective domains we have

$({\mathcal{E}}^{-1})^{\prime}(z)=\sum_{k=0}^{n}c_{k}\frac{z^{k}}{k!}+O(z^{n+1})\,,\quad
c_{k}:=\left.\frac{\mathrm{d}^{k+1}({\mathcal{E}}^{-1})}{\mathrm{d}z^{k+1}}\right|_{z=0}.$

From these formulas we obtain the estimates

$\displaystyle g(T)$
$\displaystyle=\int_{0}^{\beta{\mathcal{E}}(1)}\ln\left(1-e^{-y}\right)({\mathcal{E}}^{-1})^{\prime}(Ty)\,\mathrm{d}y$
$\displaystyle=\sum_{k=0}^{n}c_{k}T^{k}\int_{0}^{\beta{\mathcal{E}}(1)}\frac{y^{k}}{k!}\,\ln(1-e^{-y})\,\mathrm{d}y+O(T^{n+1}).$

Moreover, since

$\int_{\beta{\mathcal{E}}(1)}^{\infty}\frac{y^{k}}{k!}\,\ln(1-e^{-y})\,\mathrm{d}y=O(\beta^{k}e^{-\beta{\mathcal{E}}(1)}).$

the upper limit in all of the integrals above can be pushed to infinity at the
cost of an exponentially small term. We thus have

$\displaystyle g(T)$
$\displaystyle=\sum_{k=0}^{n}c_{k}T^{k}\int_{0}^{\infty}\frac{y^{k}}{k!}\,\ln(1-e^{-y})\,\mathrm{d}y+O(T^{n+1})$
$\displaystyle=-\sum_{k=0}^{n}c_{k}\zeta(k+2)T^{k}+O(T^{n+1})\,,$

which is equivalent to the infinite asymptotic expansion

$g(T)\sim-\sum_{k=0}^{\infty}c_{k}\zeta(k+2)T^{k}.$ (35)

By Eq. (33), the latter formula yields the following asymptotic series for the
free energy per spin of the FI and HS-B chains with $\gamma>0$:

$f_{\pm}(T)-f_{\pm}(0)\sim-\sum_{k=0}^{\infty}(m^{k+1}-1)c_{k}\zeta(k+2)T^{k+2}.$
(36)

The coefficients $c_{k}$ can be easily computed in both cases, with the result

$c_{k}=\begin{cases}(-2)^{k}(2k-1)!!\,\gamma^{-(2k+1)}&(\text{FI})\\\\[6.0pt]
(-1)^{k}(2k-1)!!\,(\gamma+1)^{-(2k+1)}&(\text{HS-B}),\end{cases}$

where $(-1)!!:=1$. In particular, from Eqs. (10)-(13) it easily follows that
the first term in the asymptotic series (36) coincides with Eq. (24).

The above argument must be slightly modified to deal with the HS chain and the
HS-B chain with $\gamma=0$, since in both cases
$({\mathcal{E}}^{-1})^{\prime}$ becomes infinite at the right endpoint
${\mathcal{E}}(1)$. For instance, for the HS chain we have

$({\mathcal{E}}^{-1})^{\prime}(z)=(1-4z)^{-1/2},\quad 0\leqslant z\leqslant
1/4;$

note that in this case ${\mathcal{E}}^{-1}$ is the inverse of
${\mathcal{E}}:[0,1/2]\to[0,1/4]$, since ${\mathcal{E}}(x)$ is increasing over
$[0,1/2]$ and symmetric about the half-point $x=1/2$. Using Eq. (27) and
performing the change of variable $\beta{\mathcal{E}}{x}=y$ we again arrive at
Eq. (33), where $g(T)$ is now given by

$g(T)=2\int_{0}^{\beta/4}\ln\left(1-e^{-y}\right)(1-4Ty)^{-1/2}\,\mathrm{d}y.$
(37)

To deal with the divergence of the last term at the upper limit of the
integral, we first note that

$(1-4z)^{-1/2}=\frac{(-1)^{n+1}}{2^{n+1}(2n+1)!!}\,\frac{\mathrm{d}^{n+1}}{\mathrm{d}z^{n+1}}(1-4z)^{\frac{2n+1}{2}}\,.$

We next define

$h(z)=\frac{(-1)^{n+1}}{2^{n+1}(2n+1)!!}\left[(1-4z)^{\frac{2n+1}{2}}-\sum_{k=0}^{2n+1}c_{k}\frac{z^{k}}{k!}\right],$

with

$c_{k}=\left.\frac{\mathrm{d}^{k}}{\mathrm{d}z^{k}}\right|_{z=0}(1-4z)^{\frac{2n+1}{2}},$

so that $h(z)=O(z^{2n+2})$. Differentiating $n+1$ times the expression for
$h(z)$ we arrive at the identity

$\displaystyle(1-4z)^{-1/2}$
$\displaystyle=h^{(n+1)}(z)+\frac{(-1)^{n+1}}{2^{n+1}(2n+1)!!}\frac{\mathrm{d}^{n+1}}{\mathrm{d}z^{n+1}}\sum_{k=n+1}^{2n+1}c_{k}\frac{z^{k}}{k!}$
$\displaystyle=h^{(n+1)}(z)+\frac{(-1)^{n+1}}{2^{n+1}(2n+1)!!}\sum_{j=0}^{n}c_{n+j+1}\frac{z^{j}}{j!}\,.$

Taking into account that

$c_{n+j+1}=(2n+1)!!(2j-1)!!(-1)^{n+1}2^{n+j+1}$

we finally obtain

$(1-4z)^{-1/2}=h^{(n+1)}(z)+\sum_{j=0}^{n}2^{j}(2j-1)!!\,\frac{z^{j}}{j!}.$

Substituting this expression into Eq. (37) yields

$\displaystyle g(T)=$
$\displaystyle\sum_{j=0}^{n}2^{j+1}(2j-1)!!\,T^{j}\int_{0}^{\beta/4}\frac{y^{j}}{j!}\ln(1-e^{-y})\mathrm{d}y$
$\displaystyle+2\beta^{n+1}\int_{0}^{\beta/4}\ln(1-e^{-y})\frac{\mathrm{d}^{n+1}}{\mathrm{d}y^{n+1}}h(Ty)\mathrm{d}y.$

The asymptotic expansion of the first term is straightforward:

$\displaystyle\sum_{j=0}^{n}2^{j+1}(2j-1)!!\,T^{j}\int_{0}^{\beta/4}\frac{y^{j}}{j!}\ln(1-e^{-y})\mathrm{d}y$
$\displaystyle=\sum_{j=0}^{n}2^{j+1}(2j-1)!!\,T^{j}\int_{0}^{\infty}\frac{y^{j}}{j!}\ln(1-e^{-y})\mathrm{d}y+O(\beta^{n}e^{-\beta/4})$
$\displaystyle=-\sum_{j=0}^{n}2^{j+1}(2j-1)!!\,\zeta(j+2)T^{j}+O(\beta^{n}e^{-\beta/4})\,.$
(38)

We claim that the second term in the previous equation for $g(T)$ is
$O(T^{n+1})$. This is easily proved by integrating by parts $n+1$ times:

$\displaystyle\beta^{n+1}\int_{0}^{\beta/4}\ln(1-e^{-y})\frac{\mathrm{d}^{n+1}}{\mathrm{d}y^{n+1}}h(Ty)\,\mathrm{d}y$
$\displaystyle=\sum_{k=0}^{n}\left.(-1)^{n-k}\beta^{n+1-k}\varphi^{(n-k)}(y)h^{(k)}(Ty)\right|^{\beta/4}_{0}$
$\displaystyle\hphantom{={}}+(-\beta)^{n+1}\int_{0}^{\beta/4}\varphi^{(n+1)}(y)h(Ty)\,\mathrm{d}y,$

with $\varphi(y):=\ln(1-e^{-y})$. The boundary terms vanish at $y=0$, since
$h(z)=O(z^{2n+2})$. On the other hand,

$\beta^{n-k+1}\varphi^{(n-k)}(\beta/4)h^{(k)}(1/4)=O(\beta^{n-k+1}e^{-\beta/4})$

for $k=0,\dots,n$, as the first $n$ derivatives of $h$ are bounded. Thus the
boundary terms are $O(\beta^{n+1}e^{-\beta/4})$, while

$\beta^{n+1}\int_{0}^{\beta/4}\varphi^{(n+1)}(y)h(Ty)\,\mathrm{d}y=O(T^{n+1}),$

since $h(Ty)=O((Ty)^{2n+2})$ and the integral
$\int_{0}^{\infty}y^{2n+2}\varphi^{(n+1)}(y)\mathrm{d}y$ is convergent.
Putting all of the above together we obtain the asymptotic series

$g(T)\sim-\sum_{k=0}^{\infty}2^{k+1}(2k-1)!!\,\zeta(k+2)T^{k},$

from which it follows that

$f_{\pm}(T)-f_{\pm}(0)\sim-\sum_{k=0}^{\infty}(m^{k+1}-1)2^{k+1}(2k-1)!!\,\zeta(k+2)T^{k+2}\,.$
(39)

In particular, truncating the series after the first term and using Eq. (10)
we obtain Eq. (28).

Consider, finally, the HS-B chain with $\gamma=0$, for which
${\mathcal{E}}(x)=x(1-x/2)$. From Eqs. (22) and (27) it easily follows that

$f_{\pm}(T)-f_{\pm}(0)=2\left(f_{\text{HS},\pm}(T/2)-f_{\text{HS},\pm}(0)\right)\,,$

whence we obtain the asymptotic series

$f_{\pm}(T)-f_{\pm}(0)\sim-\sum_{k=0}^{\infty}(m^{k+1}-1)(2k-1)!!\,\zeta(k+2)T^{k+2}\,.$

Again, the first term in this series is easily seen to yield Eq. (24).

## References

## References

  * Haldane [1988] F. D. M. Haldane, Phys. Rev. Lett. 60, 635 (1988).
  * Shastry [1988] B. S. Shastry, Phys. Rev. Lett. 60, 639 (1988).
  * Note [1] When $m=2$ we take $\gamma=1$ in the equation for $S_{i}$ to agree with the standard definition of the Pauli matrices.
  * Polychronakos [1993] A. P. Polychronakos, Phys. Rev. Lett. 70, 2329 (1993).
  * Frahm [1993] H. Frahm, J. Phys. A: Math. Gen. 26, L473 (1993).
  * Polychronakos [1994] A. P. Polychronakos, Nucl. Phys. B 419, 553 (1994).
  * Frahm and Inozemtsev [1994] H. Frahm and V. I. Inozemtsev, J. Phys. A: Math. Gen. 27, L801 (1994).
  * Note [2] In what follows sums over Latin indices will implicitly range form $1$ to $N$, unless otherwise stated.
  * Bernard  _et al._ [1995] D. Bernard, V. Pasquier, and D. Serban, Europhys. Lett. 30, 301 (1995).
  * Enciso  _et al._ [2005] A. Enciso, F. Finkel, A. González-López, and M. A. Rodríguez, Nucl. Phys. B 707, 553 (2005).
  * Laughlin [1983] R. B. Laughlin, Phys. Rev. Lett. 50, 1395 (1983).
  * Haldane [1991a] F. D. M. Haldane, Phys. Rev. Lett. 67, 937 (1991a).
  * Azuma and Iso [1994] H. Azuma and S. Iso, Phys. Lett. B 331, 107 (1994).
  * Haldane  _et al._ [1992] F. D. M. Haldane, Z. N. C. Ha, J. C. Talstra, D. Bernard, and V. Pasquier, Phys. Rev. Lett. 69, 2021 (1992).
  * Haldane [1994] F. D. M. Haldane, in  _Correlation Effects in Low-dimensional Electron Systems_ , Springer Series in Solid-state Sciences, Vol. 118, edited by A. Okiji and N. Kawakami (1994) pp. 3–20.
  * Kawakami [1992] N. Kawakami, Phys. Rev. B 46, 1005 (1992).
  * Bernard  _et al._ [1993] D. Bernard, M. Gaudin, F. D. M. Haldane, and V. Pasquier, J. Phys. A: Math. Gen. 26, 5219 (1993).
  * Ha and Haldane [1993] Z. N. C. Ha and F. D. M. Haldane, Phys. Rev. B 47, 12459 (1993).
  * Finkel and González-López [2015] F. Finkel and A. González-López, Ann. Phys.-New York 361, 520 (2015).
  * Cirac and Sierra [2010] J. I. Cirac and G. Sierra, Phys. Rev. B 81, 104431(4) (2010).
  * Nielsen  _et al._ [2011] A. E. B. Nielsen, J. I. Cirac, and G. Sierra, J. Stat. Mech.-Theory E. 2011, P11014(39) (2011).
  * Sutherland [1971] B. Sutherland, J. Math. Phys. 12, 246 (1971).
  * Sutherland [1972] B. Sutherland, Phys. Rev. A 5, 1372 (1972).
  * Ha and Haldane [1992] Z. N. C. Ha and F. D. M. Haldane, Phys. Rev. B 46, 9359 (1992).
  * Calogero [1971] F. Calogero, J. Math. Phys. 12, 419 (1971).
  * Minahan and Polychronakos [1993] J. A. Minahan and A. P. Polychronakos, Phys. Lett. B 302, 265 (1993).
  * Inozemtsev [1996] V. I. Inozemtsev, Phys. Scr. 53, 516 (1996).
  * Yamamoto [1995] T. Yamamoto, Phys. Lett. A 208, 293 (1995).
  * Finkel and González-López [2005] F. Finkel and A. González-López, Phys. Rev. B 72, 174411(6) (2005).
  * Barba  _et al._ [2008] J. C. Barba, F. Finkel, A. González-López, and M. A. Rodríguez, Europhys. Lett. 83, 27005(6) (2008).
  * Barba  _et al._ [2010] J. C. Barba, F. Finkel, A. González-López, and M. A. Rodríguez, Nucl. Phys. B 839, 499 (2010).
  * Kirillov  _et al._ [1997] A. N. Kirillov, A. Kuniba, and T. Nakanishi, Commun. Math. Phys. 185, 441 (1997).
  * Basu-Mallick  _et al._ [2010] B. Basu-Mallick, N. Bondyopadhaya, and K. Hikami, Symmetry Integr. Geom. 6, 091(13) (2010).
  * Carrasco  _et al._ [2020] J. Carrasco, F. Finkel, A. González-López, and M. A. Rodríguez (2020), arXiv:2010.08307 [math-ph], to appear in Adv. Theor. Math. Phys.
  * Haldane [1991b] F. D. M. Haldane, Phys. Rev. Lett. 66, 1529 (1991b).
  * Sutherland and Shastry [1993] B. Sutherland and B. S. Shastry, Phys. Rev. Lett. 71, 5 (1993).
  * Enciso  _et al._ [2012] A. Enciso, F. Finkel, and A. González-López, Ann. Phys.-New York 327, 2627 (2012).
  * Finkel  _et al._ [2018] F. Finkel, A. González-López, I. León, and M. A. Rodríguez, J. Stat. Mech.-Theory E. 2018, 043101(47) (2018).
  * Perron [1907] O. Perron, Math. Ann. 64, 1 (1907).
  * Gantmacher and Krein [2000] F. Gantmacher and M. Krein,  _Oscillation Matrices and Kernels and Small Vibrations of Mechanical Systems_ , revised ed. (AMS Chelsea Publishing, Providence, RI, 2000).
  * Finkel and González-López [2022] F. Finkel and A. González-López (2022), arXiv:2206.02651v1 [cond-mat.stat-mech], to appear in J. Stat. Mech.-Theory E.
  * Note [3] In the truly supersymmetric case the transfer matrix has always a zero eigenvalue, which is doubly degenerate for $m=n=2$. This makes it straightforward to diagonalize the latter matrix when $1\leqslant m,n\leqslant 2$.
  * Blöte  _et al._ [1986] H. W. J. Blöte, J. L. Cardy, and M. P. Nightingale, Phys. Rev. Lett. 56, 742 (1986).
  * Affleck [1986] I. Affleck, Phys. Rev. Lett. 56, 746 (1986).
  * Hikami and Basu-Mallick [2000] K. Hikami and B. Basu-Mallick, Nucl. Phys. B 566, 511 (2000).
  * Schoutens [1994] K. Schoutens, Phys. Lett. B 331, 335 (1994).
  * Bouwknegt and Schoutens [1996] P. Bouwknegt and K. Schoutens, Nucl. Phys. B 482, 345 (1996).
  * Basu-Mallick  _et al._ [2008] B. Basu-Mallick, N. Bondyopadhaya, and D. Sen, Nucl. Phys. B 795, 596 (2008).
  * Note [4] A matrix (or, in particular, a vector) is said to be positive if all its entries entries are positive.
  * Note [5] For the sake of conciseness, we have omitted the dependence of $a$ on $x$.
  * Olver  _et al._ [2010] F. W. J. Olver, D. W. Lozier, R. F. Boisvert, and C. W. Clark, eds.,  _NIST Handbook of Mathematical Functions_ (Cambridge University Press, 2010).
  * Note [6] In fact, it is shown in the appendix that for the PF chain the $O(T^{3})$ term is actually $O(T^{2}e^{-\beta/m})$.
  * Note [7] We use the notation $\lfloor x\rfloor$ to denote the integer part of the real number $x$.
  * Hikami [1995] K. Hikami, J. Phys. Soc. Jpn. 64, 1047 (1995).


# Oscillatory Edge Modes in Two Dimensional Spin-Torque Oscillator Arrays

Shivam Kamboj  Department of Physics and Astronomy, San José State University,
San José, California, 95192, USA  Department of Physics, University of
Califonia, Merced, California, 95343, USA  Rembert A. Duine  Institute for
Theoretical Physics and Center for Extreme Matter and Emergent Phenomena,
Utrecht University, Leuvenlaan 4, 3584 CE Utrecht, The Netherlands  Department
of Applied Physics, Eindhoven University of Technology, P.O. Box 513, 5600 MB
Eindhoven, The Netherlands  Benedetta Flebus  Department of Physics, Boston
College, 140 Commonwealth Avenue, Chestnut Hill, Massachusetts 02467, USA
Hilary M. Hurst  Department of Physics and Astronomy, San José State
University, San José, California, 95192, USA

###### Abstract

Spin torque oscillators (STOs) are dissipative magnetic systems that provide a
natural platform for exploring non-Hermitian phenomena. We theoretically study
a two-dimensional (2d) array of STOs and show that its dynamics can be mapped
to a 2d, non-Hermitian Su-Schrieffer-Heeger (SSH) model. We calculate the
energy spectrum and identify the one-dimensional (1d) edge states of our
model, corresponding to auto-oscillation of STOs on the boundary of the system
while the bulk oscillators do not activate. We show that tuning the Gilbert
damping, injected spin current, and coupling between STOs allows for exploring
the edge state properties under different parameter regimes. Furthermore, this
system admits 1d edge states with non-uniform probability density, and we
explore their properties in systems of different sizes. Additional symmetry
analysis indicates that these states are not topologically protected but are
nevertheless confined to the edge of the system, as the bulk is protected by
$\mathcal{P}\mathcal{T}$-symmetry. These results indicate that 2d arrays of
STOs may be useful to explore novel edge state behavior in dissipative
systems.

##  I Introduction

Topology and its connection to condensed matter systems has been the subject
of intense research for nearly fifty years, since the discovery of quantized
Hall resistance and its topological origin [1, 2]. Topology is now understood
as a critical underlying feature of many materials that can affect global
transport properties, resulting in, e.g., the quantum anomalous and spin Hall
effects [3, 4] and leading to entirely new classes of topological materials
[5, 6, 7]. The effects of topology in systems with dissipation, i.e. non-
Hermitian systems, can be markedly different from their Hermitian counterparts
[8, 9, 10, 11, 12, 13, 14]. Non-Hermitian systems can exhibit exceptional
points [15, 16, 17, 18], i.e. the coalescence of two or more eigenvectors, as
well as the non-Hermitian skin effect, a phenomenon where bulk eigenstates
localize on the edge of the system [19, 20, 21]. The edge states in non-
Hermitian systems can also exist despite the breakdown of the bulk-boundary
correspondence [9, 10] and exhibit lasing behavior [22, 23]. While several
experimental realizations of non-Hermitian phenomena have been observed in
photonic [22], acoustic [24], and electronic circuits [25, 26], their
experimental exploration in magnonic systems is still in its infancy. Magnonic
systems are, however, a natural platform in which to realize non-Hermitian
physics because they are always coupled to a surrounding environment and
exhibit lossy dynamics [27, 28].

Spin torque oscillators (STOs) have recently emerged as a promising platform
for harboring non-Hermitian phenomena [29, 27, 30]. These magnetic nanopillars
are nanometer sized devices that conduct spin currents via spin transfer
torque [31]. STOs are dissipative systems because, like all magnetic systems,
they are subject to ubiquitous spin non-conserving interactions parameterized
by Gilbert damping [32]. The magnetic dynamics of 1d STO arrays was
successfully mapped to a non-Hermitian Su-Schrieffer-Heeger (SSH) model with
topologically protected lasing edge states in Ref. [29]. In a subsequent
numerical study these edge states were shown to be robust in the presence of
additional terms such as dipolar interactions and nonlinear STO behavior [33].
Here, we examine whether these lasing edge states can be realized in higher-
dimensional arrays of STOs.

The 1d SSH model and it’s non-Hermitian variants are widely used in condensed
matter physics to study topological systems [34, 10]; 2d non-Hermitian SSH
models have been studied in a variety of contexts, where a number of
interesting properties such as in-gap topological states and non-trivial bulk-
band topology have been found [35, 36, 37, 25, 26]. However, most of these
models lack a clear experimental implementation, with topological circuits
being a notable exception [25, 26]. The sources of non-Hermitian terms (i.e.
energy non-conserving terms) can be difficult to quantify and characterize in
experimental platforms. Here, we propose STO arrays as a platform for
experimental realization of a non-Hermitian 2d SSH model. Specifically, we
focus on a geometry consisting of several 1d chains that are weakly coupled to
form a 2d STO array. By introducing this new type of ‘vertical coupling’, we
derive a non-Hermitian 2d SSH model that exhibits 1d lasing edge states.

The manuscript is organized as follows: In Section II we map the linearized
Landau-Lifshitz-Gilbert (LLG) equation for the magnetization dynamics into a
non-Hermitian tight-binding Hamiltonian. In Sec. III we discuss the properties
of our model including the oscillatory edge states and the symmetry properties
of the Bloch Hamiltonian. Finally, we conclude and discuss directions for
future work in Sec. IV.

##  II Model

We consider a 2d array of $M\times 2N$ STOs which could be fabricated from
individual nanopillars [29] or a multilayer structure [33]. Here, $M$ is the
total number of rows in the array and $N$ indicates the number of unit cells
per row, where there are two STOs per unit cell. A single STO consists of a
layer of fixed magnetic polarization and a ‘free’ magnetic layer without fixed
polarization, separated by a thin metallic spacer. By injecting spin current
into the free layer, the fixed layer is driven to precess about its
equilibrium direction, which is set by an external applied magnetic field. The
dynamics of an isolated STO subjected to a magnetic field
$\mathbf{H}_{0}=H_{0}\hat{\mathbf{z}}$ and spin current
$\mathbf{J}^{S}_{\eta}=J^{S}_{\eta}\hat{\mathbf{z}}$ are described by the LLG
equation for the magnetization vector $\mathbf{m}_{\eta,ij}$ [31]. Here, the
index $\eta=A,B$ denotes the $A$ and $B$ sublattices and the indices $i,j$
label the sites using the (row, column) convention. The resulting LLG equation
is

$\displaystyle\dot{\mathbf{m}}_{\eta,ij}\left|{}_{0}\right.$
$\displaystyle=\omega_{\eta,ij}\hat{\mathbf{z}}\times\mathbf{m}_{\eta,ij}+\alpha_{\eta,ij}\mathbf{m}_{\eta,ij}\times\dot{\mathbf{m}}_{\eta,ij}$
$\displaystyle+J^{S}_{\eta}\mathbf{m}_{\eta,ij}\times(\mathbf{m}_{\eta,ij}\times\hat{\mathbf{z}}).$
(1)

The STO ferromagnetic resonance frequency is given by
$\omega_{\eta,ij}=\gamma_{\eta,ij}(H_{0}-4\pi M_{\eta,ij})$ where
$\gamma_{\eta,ij}$ is the gyromagnetic ratio and $M_{\eta,ij}$ is the
saturation magnetization; $\alpha_{\eta,ij}\ll 1$ is the dimensionless Gilbert
damping parameter. We assume the resonance frequency and the Gilbert damping
to the be same for all STOs in the array and drop the subscripts going
forward, i.e. $\omega_{\eta,ij}\rightarrow\omega$ and
$\alpha_{\eta,ij}\rightarrow\alpha$. The third term in Eq. (1) is the
Slonczewski-Berger spin-transfer torque [38, 39]. The injected spin current
$J^{S}_{\eta}$ is assumed to the be same for all sites in a given sublattice.

Interactions between adjacent STOs are mediated by a reactive Ruderman-
Kittel-Kasuya-Yosida (RKKY)-type magnetic exchange coupling, which is given by
[29]

$\displaystyle\dot{\mathbf{m}}_{A,ij}\left|{}_{\rm coup}\right.$
$\displaystyle=-\mathbf{m}_{A,ij}\times(J\mathbf{m}_{B,ij}+\tilde{J}\mathbf{m}_{B,ij-1})$
$\displaystyle-\mathbf{m}_{A,ij}\times
J_{2}(\mathbf{m}_{\eta^{\prime},i+1j}+\mathbf{m}_{\eta^{\prime},i-1j})$ (2)

$\displaystyle\dot{\mathbf{m}}_{B,ij}\left|{}_{\rm
coup}\right.=-\mathbf{m}_{B,ij}\times(J\mathbf{m}_{A,ij}+\tilde{J}\mathbf{m}_{A,ij+1})$
$\displaystyle-\mathbf{m}_{B,ij}\times
J_{2}(\mathbf{m}_{\eta^{\prime\prime},i+1j}+\mathbf{m}_{\eta^{\prime\prime},i-1j})$
(3)

The coupling strengths $J,\tilde{J}$ are the intracell and intercell coupling
for the 1d unit cell, which contains two STOs, and $J,\tilde{J},J_{2}>0$
indicate ferromagnetic coupling. We consider the geometry depicted in Fig. 1,
for which $\eta^{\prime}=A,\eta^{\prime\prime}=B$.

Figure 1: Schematic of the effective tight-binding model of a 2d STO array,
showing a system of size $4\times 8$, where there are $4$ unit cells per row.
Red (dark color) sites are the A sublattice and green (light color) sites are
the B lattice.

We now derive the linearized equations of motion for the STO magnetization
$\mathbf{m}_{A(B)}$. That is, we start in the strong-field regime where the
magnetic moment is mostly aligned along the $\hat{\mathbf{z}}$ direction and
the energy scale given by $\omega$ is the largest one in the problem. Starting
from Eq. (1), we can derive an effective Hamiltonian in the following way:
First, we linearize the equations of motion about the equilibrium field
direction by writing
$\mathbf{m}_{\eta,ij}=(m^{x}_{\eta,ij},m^{y}_{\eta,ij},1)^{T}$ where
$|m^{x}_{\eta,ij}|\simeq|m^{y}_{\eta,ij}|\ll 1$, and $|m^{z}|\simeq 1$ is
assumed to be a constant. We then introduce the variable
$2m^{-}_{\eta,ij}=m^{x}_{\eta,ij}-im^{y}_{\eta,ij}$, and make the Holstein-
Primakoff approximation $m^{-}_{A(B),ij}=\langle a_{ij}(b_{ij})\rangle
e^{-i\omega t}$ where the second quantized bosonic operators $a_{ij}$,
$b_{ij}$ annihilate a magnon on the STO at site $(i,j)$ [40]. From the
Heisenberg equations of motion $\dot{a}_{ij}=i/\hbar\left[H,a_{ij}\right]$, we
derive an effective Hamiltonian corresponding to the linearized equations of
motion. In the following we set $\hbar=1$.

Figure 2: The dependence of the real (a-b) and imaginary (e-f) parts of the
energy spectrum of $H$ on $J/|\tilde{J}|$ for a system of $10\times 20$ STOs
for $\omega=1\tilde{J}$, $\alpha=0.2$, and $J_{2}=0.01\tilde{J}$ (a,e) and
$J_{2}=0.1\tilde{J}$ (b,f). Modes with $\mathrm{Im}(E)>0$ correspond to auto-
oscillation of STOs on the edge of the system, appearing for $|J|<\tilde{J}$.
These modes correspond to the flat in-gap bands at
$\mathrm{Re}(E)-\omega\approx 0$ in the real spectrum. Panels (c-d, g-h) show
the real (c-d) and imaginary (g-h) parts of the energy spectrum for a system
of $50\times 100$ STOs, with $J_{2}=0.01\tilde{J}$ (c,g) and
$J_{2}=0.1\tilde{J}$ (d,h). Apart from the different system size, all the
other parameters are same as in (a-b) and (e-f). In the larger system there is
some activation of bulk STOs in the $\mathcal{P}\mathcal{T}$-broken regime,
indicated in (g-h) by the additional states with $\mathrm{Im}(E)>0$ in the
region $\alpha\omega>|\tilde{J}-J|$.

For the 2d STO array, $H=\sum_{ij}H_{ij}$ is

$\displaystyle H_{ij}$
$\displaystyle=\omega(a^{\dagger}_{ij}a_{ij}+b^{\dagger}_{ij}b_{ij})$
$\displaystyle+i(J^{S}_{A}-\alpha\omega)a^{\dagger}_{ij}a_{ij}+i(J^{S}_{B}-\alpha\omega)b^{\dagger}_{ij}b_{ij}$
$\displaystyle-J(a^{\dagger}_{ij}b_{ij}+h.c.)-\tilde{J}(a^{\dagger}_{ij}b_{ij-1}+h.c.)$
$\displaystyle-J_{2}(a^{\dagger}_{ij}a_{i-1j}+b^{\dagger}_{ij}b_{i-1j}+h.c).$
(4)

We see that non-Hermiticity arises due to the onsite spin current injection
and Gilbert damping, resulting in onsite terms $\propto
i(J^{S}_{\eta}-\alpha\omega)$, and the degree of non-Hermiticity can be tuned
by balancing the injected spin current and Gilbert damping. Given the non-
Hermitian lattice model, we can now explore its energy spectrum and edge state
properties. In this manuscript we use the term ‘energy spectrum’ and the
symbol $E$ to denote the complex eigenvalues of $H$; Re($E$) can be thought of
as an energy while $\mathrm{Im}(E)>0$ ($<0$) is an indication of lasing
(damping).

Here we consider the effect of different coupling strengths $J_{2}$. We
confine the system to the parity-time ($\mathcal{P}\mathcal{T}$) symmetric
regime where the injected spin current is $J_{SA}=2\alpha\omega$, $J_{SB}=0$;
further symmetry analysis is performed in Sec. III.3. This is a 2d extension
of the 1d $\mathcal{P}\mathcal{T}$-symmetric, non-Hermitian SSH model analyzed
in Ref. [29]. We note that other versions of the non-Hermitian 2d SSH model
studied in the literature have alternating A and B sublattice sites in the
vertical direction [35, 36, 37], resulting in a four site unit cell. Here,
however, we introduce the coupling such that each column is all A or B
sublattice sites, thus reducing the unit cell to two STOs. Our reasoning for
this geometry is that it would be easier to inject spin current into an entire
column of STOs, for example using a metallic strip, rather than having to
individually address each A or B site STO in the grid.

##  III Results

Figure 3: Spatial distribution of 1d edge modes in a 2d STO array. (a-d) Color
density plots showing $|\psi|^{2}$ for four of the the flat-band modes with
$\mathrm{Im}(E)>0$. The states are localized on the left edge, with
$|\psi|^{2}=0$ everywhere in the bulk. The insets show the 1d distribution
along the edge. (e) 1d plot showing the spatial distribution of all 10 lasing
edge states and the corresponding real part of the energy,
$\mathrm{Re}(E)-\omega$. The edge states are not uniform across the system due
to the coupling $J_{2}$ between rows. States appear in pairs where states with
equal magnitude and opposite sign of $\mathrm{Re}(E)-\omega$ have the same
spatial distribution (dashed and solid lines overlap).

Our goal is to understand how the vertical coupling $J_{2}$ between rows
affects the properties of the original 1d system, thus we briefly review the
results for a 1d STO chain [29, 33]. The real energy spectrum of the 1d STO
model Hamiltonian has two degenerate flat bands for $J<|\tilde{J}|$ and admits
a real line gap in $k$-space [12]. The flat bands with degenerate energy
eigenvalues are an indication of topologically protected edge states of the
Hamiltonian. The eigenstates corresponding to the flat-bands also have non-
zero imaginary eigenvalues, indicating lasing ($\mathrm{Im}(E)>0$) and damped
($\mathrm{Im}(E)<0$) states. The lasing states correspond to an auto
oscillation of STOs, which occurs only at the edge of the system for
$\alpha\omega<|J-\tilde{J}|$. In the regime $\alpha\omega>|J-\tilde{J}|$, the
bulk Hamiltonian also has complex eigenvalues which can lead to oscillation of
the bulk STOs; this is the so-called $\mathcal{P}\mathcal{T}$-broken regime of
the model where the Hamiltonian respects $\mathcal{P}\mathcal{T}$-symmetry but
its eigenstates do not [27]. Symmetry analysis of the Bloch Hamiltonian for
the 1d model confirms the auto-oscillation of the edge STOs to be a zero-
dimensional ‘edge state’ with topological protection.

In the 2d case, the non-Hermitian Hamiltonian Eq. (4) also has complex
eigenvalues. The energy spectrum of the 2d model looks similar to the 1d model
in the case of weak vertical coupling $J_{2}/\tilde{J}\lesssim 0.01$, but the
degeneracy in the flat bands breaks immediately even with infinitesimal
vertical coupling. In Fig. 2, we show the energy spectrum of the Hamiltonian
considering different vertical coupling strengths $J_{2}/\tilde{J}=0.01,0.1$
for two different system sizes. We found that as $J_{2}$ increases with
respect to $\tilde{J}$, the energy separation between flat bands increases and
they hybridize with the bulk states. We simulated results for a 2d array 20
sites wide and 10 sites in the vertical direction, as well as a larger system
100 sites wide and 50 sites in the vertical direction. In the upcoming
subsections, we respectively discuss the numerical results for the edge
states, analyze the Hamiltonian in momentum space, and present a symmetry
analysis of the model.

###  III.1 Edge States

Using exact diagonalization, we find that the real-space Hamiltonian Eq. (4)
exhibits one-dimensional ‘lasing’ edge states where $\mathrm{Im}(E)>0$. In
Fig. 3 we show the spatial distribution of edge states for a system of
$10\times 20$ STOs. Since there are 10 rows, the system exhibits 10 lasing
edge modes, all with $\mathrm{Im}(E)=0.09$. We have confirmed that the
corresponding damped edge modes with $\mathrm{Im}(E)=-0.09$ occur on the
opposite edge of the system, as expected (not pictured here). We set the
intercell coupling to $\tilde{J}=1$ and consider the regime where intracell
coupling $J/\tilde{J}=0.2$ and vertical coupling $J_{2}/\tilde{J}=0.1$. The
physical manifestation of these edge states is an auto-oscillation of STOs
that is non-uniform at the edge of the sample. Thus, the STOs will exhibit
spatially varying microwave emission which can be tuned based on the magnon
population in each edge mode.

Figure 4: Example of the spatial distribution of 10 unique lasing edge states
(out of 50) and the corresponding real part of the energy,
$\mathrm{Re}(E)-\omega\times 10^{-1}$ for a system of $50\times 100$ sites.
The larger system shows greater heterogeneity in the spatial distribution of
edge modes. The top row of the legend corresponds to the left column and the
bottom row corresponds to the right column; values increase downward in each
column. For example $\mathrm{Re}(E)-\omega=1.0$ corresponds to the bottom-left
panel.

 We also diagonalized a larger STO array to see how the results vary if we
scale up the system size. The Hamiltonian matrix has a block banded structure,
called a block Toeplitz-tridiagonal (TT) matrix. For a larger system i.e.
50x100 sites, the matrix becomes 5000 x 5000 for which the diagonalization is
computationally expensive for different values of $J_{2}$. To diagonalize this
large sparse matrix efficiently, we used the method given by Ref. [41]. The
energy spectrum for the larger system - as shown in Fig. 2 (c-d), (g-h) -
doesn’t show any deviation apart from very small bulk oscillations in the
region where $\alpha\omega>|J-\tilde{J}|$. This proves the 1d edge states can
exist in a significantly larger system. As shown in Fig. 4, the edge states in
the larger system also exhibit a non-uniform spatial distribution, with
increased oscillations and beating behavior visible in some modes. Here, we
show the first 10 out of 50 total edge states ranked by increasing Re($E$).

###  III.2 Bloch Hamiltonian

Here we consider the $\mathcal{P}\mathcal{T}$-symmetric regime where
$J_{SA}=2\alpha\omega$ and $J_{SB}=0$, leading to balanced gain (loss) terms
$\pm i\alpha\omega$ on the A(B) sites. To analyze the Hamiltonian in momentum
space, we consider periodic boundary conditions and Fourier transform Eq. (4),
which is written as

$\hat{H}=\sum_{\mathbf{k}}\begin{pmatrix}a^{\dagger}_{\mathbf{k}}&b^{\dagger}_{\mathbf{k}}\end{pmatrix}\left(H_{\mathbf{k}}-\omega\right)\begin{pmatrix}a_{\mathbf{k}}\\\
b_{\mathbf{k}}\end{pmatrix}$ (5)

where

$H_{\mathbf{k}}=\left(\begin{array}[]{cc}i\alpha\omega-2J_{2}\cos{k_{y}}&-J-\tilde{J}e^{ik_{x}}\\\
-J-\tilde{J}e^{-ik_{x}}&-i\alpha\omega-2J_{2}\cos{k_{y}}\end{array}\right)$
(6)

The resonant frequency simply provides an overall shift of the energy
spectrum, therefore we redefine $\omega$ as the zero energy point.

The Bloch Hamiltonian can be written

$H_{\mathbf{k}}=d_{0}(k_{y})\mathbb{1}+\mathbf{d}(k_{x})\cdot\bm{\sigma}$ (7)

where $\mathbb{1}$ is the $2\times 2$ identity matrix and
$\bm{\sigma}=(\sigma^{x},\sigma^{y},\sigma^{z})$ is the vector of Pauli
matrices. We define the functions $d_{0}(k_{y})=-2J_{2}\cos(k_{y})$ and

$\mathbf{d}(k_{x})=\begin{pmatrix}-J-\tilde{J}\cos(k_{x})\\\
\tilde{J}\sin(k_{x})\\\ i\alpha\omega\end{pmatrix}.$ (8)

We find the two-band energy spectrum

$\epsilon_{\pm}(\mathbf{k})=-2J_{2}\cos(k_{y})\pm\sqrt{J^{2}+\tilde{J}^{2}+2J\tilde{J}\cos(k_{x})-\alpha^{2}\omega^{2}}.$
(9)

Here we see that the eigenvalues are real for
$\alpha\omega<|J+\tilde{J}e^{ik_{x}}|$, i.e. the system remains in the
$\mathcal{P}\mathcal{T}$-unbroken regime exhibiting real eigenvalues as long
as the Gilbert damping is relatively small. This condition is satisfied for
all $k$ if $\alpha\omega<|J-\tilde{J}|$. Furthermore, the spectrum is linear
in the vertical coupling $J_{2}$; thus the effect of coupling adjacent 1d STO
chains together is to shift the spectrum away from the resonant frequency
$\omega$. We can see this clearly for example in Fig. 2 (b), where the flat
bands of adjacent 1d chains hybridize and are vertically shifted. Furthermore,
$J_{2}$ can cause the energy gap to close, however in the regime
$J_{2}\lesssim\sqrt{(J-\tilde{J})^{2}-\alpha^{2}\omega^{2}}$ the energy gap is
open and the edge states remain well separated from the bulk.

###  III.3 Symmetry Analysis

Symmetry analysis can help determine whether the 1d edge states displayed in
Figs. 3 and 4 are topologically protected. We investigate the following
symmetries of the Bloch Hamiltonian, Eq. (6): chiral symmetry, chiral-
inversion symmetry, sublattice symmetry, and parity-time
($\mathcal{P}\mathcal{T}$) symmetry. Systems obeying chiral, chiral-inversion,
or sublattice symmetry can exhibit topologically protected edge modes [11, 42,
12], therefore it is important to check whether these symmetries are preserved
for our model. $\mathcal{P}\mathcal{T}$-symmetry ensures there is always a
regime in which the Hamiltonian has real eigenvalues [43]. We note that for a
Hermitian system, chiral and sublattice symmetries are equivalent, however for
a non-Hermitian Hamiltonian this is no longer the case and some care must be
taken. Here we use the symmetry naming conventions from Ref. [12].

To have chiral symmetry (CS), the Hamiltonian must satisfy the condition
$\sigma_{z}H_{\mathbf{k}}^{\dagger}\sigma_{z}=-H_{\mathbf{k}}$. We find that
the vertical coupling term $J_{2}\cos(k_{y})$ breaks chiral symmetry in
general, however for the special values $k_{y}=\pm\pi/2$ chiral symmetry is
preserved. To have chiral-inversion (CI) symmetry, the Hamiltonian must
satisfy the condition $\sigma_{y}H_{\mathbf{k}}\sigma_{y}=-H_{-\mathbf{k}}$.
Like the case of CS, CI is in general broken by $J_{2}$ and only preserved for
$k_{y}=\pm\pi/2$. For sublattice symmetry, the Hamiltonian must satisfy the
condition $\sigma_{z}H_{\mathbf{k}}\sigma_{z}=-H_{\mathbf{k}}$. This model
does not posses sublattice symmetry for any parameter regime due to the non-
Hermitian terms, as is the case in 1d [10, 29].

To have $\mathcal{P}\mathcal{T}$-symmetry, the Hamiltonian must satisfy the
condition $\sigma_{x}H_{\mathbf{k}}^{*}\sigma_{x}=H_{\mathbf{k}}$. This
condition is satisfied for $J_{SA}=2\alpha\omega$ and $J_{SB}=0$, thus as with
the 1d case the system can be tuned to the $\mathcal{P}\mathcal{T}$-symmetric
regime by altering the injected spin current on A and B sublattice sites.
$\mathcal{P}\mathcal{T}$-symmetry alone does not guarantee topological
protection [10], and the results from the symmetry analysis indicate that the
edge states observed in this model are not topologically protected. However,
the confinement of the oscillatory modes to the edge can be understood as a
result of $\mathcal{P}\mathcal{T}$-symmetry in the bulk, which is broken
spontaneously by the edge of the system. Furthermore,
$\mathcal{P}\mathcal{T}$-symmetry guarantees that the bulk STOs do not have
any lasing modes as long as $\alpha\omega<|J-\tilde{J}|$.

##  IV Conclusion and Outlook

In this work we have examined a novel realization of a non-Hermitian 2d SSH
model which can be constructed from an array of STOs. Using exact
diagonalization and analysis of the Bloch Hamiltonian, we have shown that this
model exhibits 1d lasing edge states with a non-uniform spatial distribution.
The physical manifestation of these modes is an auto-oscillation of STOs along
one edge of the system which is spatially varying. The extension of the model
from 1d to 2d via the addition of vertical coupling between individual 1d STO
chains breaks chiral-inversion, chiral, and sublattice symmetry, indicating a
loss of topological protection for these modes. However, the vertical coupling
preserves $\mathcal{P}\mathcal{T}$-symmetry for the bulk states, thereby
guaranteeing that the bulk oscillators do not activate, even in the presence
of spin current injected into the bulk.

Here we have considered an injected spin current such that the system remains
at the $\mathcal{P}\mathcal{T}$-symmetric point. Future works could
investigate the robustness of these edge states in the presence of additional
terms such as dipolar interactions as well as dissipative coupling between
STOs. However, results from studies of the analogous 1d model indicate that if
these terms are small compared to the reactive RKKY coupling studied here,
they would not strongly affect the presence of edge states [29, 33].

Another interesting extension of this model would be to explore possible
application in devices. The dispersion relation in Eq. (9) indicates that the
1d edge states have nonzero group velocity due to the coupling $J_{2}$. Future
work can examine whether these edge states, despite being non-uniform, could
be used as a 1d channel for spin transport. For example, if there is
additional spin injected into one of the edge STOs exhibiting lasing, one
could study how it travels along the edge of the system. Such a device could
potentially provide a new way to realize a low-dissipation transport channel
for spin.

###### Acknowledgements.

H.M.H. acknowledges support of the San José State University Research,
Scholarship, and Creative Activity assigned time program. B.F. acknowledges
support of the National Science Foundation under Grant No. NSF DMR-2144086.

## References

  * Klitzing  _et al._ [1980] K. v. Klitzing, G. Dorda, and M. Pepper, New method for high-accuracy determination of the fine-structure constant based on quantized hall resistance, Physical review letters 45, 494 (1980).
  * Thouless  _et al._ [1982] D. J. Thouless, M. Kohmoto, M. P. Nightingale, and M. den Nijs, Quantized hall conductance in a two-dimensional periodic potential, Physical review letters 49, 405 (1982).
  * Maciejko  _et al._ [2011] J. Maciejko, T. L. Hughes, and S.-C. Zhang, The quantum spin hall effect, Annu. Rev. Condens. Matter Phys. 2, 31 (2011).
  * Liu  _et al._ [2016] C.-X. Liu, S.-C. Zhang, and X.-L. Qi, The quantum anomalous hall effect: theory and experiment, Annual Review of Condensed Matter Physics 7, 301 (2016).
  * Moore [2010] J. E. Moore, The birth of topological insulators, Nature 464, 194 (2010).
  * Sato and Ando [2017] M. Sato and Y. Ando, Topological superconductors: a review, Reports on Progress in Physics 80, 076501 (2017).
  * Yan and Felser [2017] B. Yan and C. Felser, Topological materials: Weyl semimetals, Annual Review of Condensed Matter Physics 8, 337 (2017).
  * Hu and Hughes [2011] Y. C. Hu and T. L. Hughes, Absence of topological insulator phases in non-hermitian $\mathcal{P}\mathcal{T}$-symmetric hamiltonians, Physical Review B 84, 153101 (2011).
  * Yao and Wang [2018] S. Yao and Z. Wang, Edge states and topological invariants of non-hermitian systems, Phys. Rev. Lett. 121, 086803 (2018).
  * Lieu [2018] S. Lieu, Topological phases in the non-hermitian su-schrieffer-heeger model, Physical Review B 97, 045106 (2018).
  * Gong  _et al._ [2018] Z. Gong, Y. Ashida, K. Kawabata, K. Takasan, S. Higashikawa, and M. Ueda, Topological phases of non-hermitian systems, Physical Review X 8, 031079 (2018).
  * Kawabata  _et al._ [2019] K. Kawabata, K. Shiozaki, M. Ueda, and M. Sato, Symmetry and topology in non-hermitian physics, Physical Review X 9, 041015 (2019).
  * Bergholtz  _et al._ [2021] E. J. Bergholtz, J. C. Budich, and F. K. Kunst, Exceptional topology of non-hermitian systems, Rev. Mod. Phys. 93, 015005 (2021).
  * Ding  _et al._ [2022] K. Ding, C. Fang, and G. Ma, Non-hermitian topology and exceptional-point geometries, Nature Reviews Physics , 1 (2022).
  * Heiss [2012] W. Heiss, The physics of exceptional points, Journal of Physics A: Mathematical and Theoretical 45, 444016 (2012).
  * Doppler  _et al._ [2016] J. Doppler, A. A. Mailybaev, J. Böhm, U. Kuhl, A. Girschik, F. Libisch, T. J. Milburn, P. Rabl, N. Moiseyev, and S. Rotter, Dynamically encircling an exceptional point for asymmetric mode switching, Nature 537, 76 (2016).
  * Liu  _et al._ [2019a] H. Liu, D. Sun, C. Zhang, M. Groesbeck, R. Mclaughlin, and Z. V. Vardeny, Observation of exceptional points in magnonic parity-time symmetry devices, Science advances 5, eaax9144 (2019a).
  * Deng  _et al._ [2023] K. Deng, X. Li, and B. Flebus, Exceptional points as signatures of dynamical magnetic phase transitions, Physical Review B 107, L100402 (2023).
  * Yokomizo and Murakami [2021] K. Yokomizo and S. Murakami, Scaling rule for the critical non-hermitian skin effect, Physical Review B 104, 165117 (2021).
  * Deng and Flebus [2022] K. Deng and B. Flebus, Non-hermitian skin effect in magnetic systems, Physical Review B 105, L180406 (2022).
  * Zhang  _et al._ [2022] K. Zhang, Z. Yang, and C. Fang, Universal non-hermitian skin effect in two and higher dimensions, Nature communications 13, 2496 (2022).
  * Feng  _et al._ [2017] L. Feng, R. El-Ganainy, and L. Ge, Non-hermitian photonics based on parity–time symmetry, Nature Photonics 11, 752 (2017).
  * Ota  _et al._ [2020] Y. Ota, K. Takata, T. Ozawa, A. Amo, Z. Jia, B. Kante, M. Notomi, Y. Arakawa, and S. Iwamoto, Active topological photonics, Nanophotonics 9, 547 (2020).
  * Zhu  _et al._ [2014] X. Zhu, H. Ramezani, C. Shi, J. Zhu, and X. Zhang, $\mathcal{P}\mathcal{T}$-symmetric acoustics, Phys. Rev. X 4, 031042 (2014).
  * Helbig  _et al._ [2020] T. Helbig, T. Hofmann, S. Imhof, M. Abdelghany, T. Kiessling, L. Molenkamp, C. Lee, A. Szameit, M. Greiter, and R. Thomale, Generalized bulk–boundary correspondence in non-hermitian topolectrical circuits, Nature Physics 16, 747 (2020).
  * Kotwal  _et al._ [2021] T. Kotwal, F. Moseley, A. Stegmaier, S. Imhof, H. Brand, T. Kießling, R. Thomale, H. Ronellenfitsch, and J. Dunkel, Active topolectrical circuits, Proceedings of the National Academy of Sciences 118, e2106411118 (2021).
  * Hurst and Flebus [2022] H. M. Hurst and B. Flebus, Non-hermitian physics in magnetic systems, Journal of Applied Physics 132, 220902 (2022).
  * Yuan  _et al._ [2022] H. Yuan, Y. Cao, A. Kamra, R. A. Duine, and P. Yan, Quantum magnonics: when magnon spintronics meets quantum information science, Physics Reports 965, 1 (2022).
  * Flebus  _et al._ [2020] B. Flebus, R. A. Duine, and H. M. Hurst, Non-hermitian topology of one-dimensional spin-torque oscillator arrays, Phys. Rev. B 102, 180408 (2020).
  * Wittrock  _et al._ [2023] S. Wittrock, S. Perna, R. Lebrun, K. Ho, R. Dutra, R. Ferreira, P. Bortolotti, C. Serpico, and V. Cros, Non-hermiticity in spintronics: oscillation death in coupled spintronic nano-oscillators through emerging exceptional points (2023), arXiv:2108.04804 [cond-mat.mes-hall] .
  * Slavin and Tiberkevich [2009] A. Slavin and V. Tiberkevich, Nonlinear auto-oscillator theory of microwave generation by spin-polarized current, IEEE Transactions on Magnetics 45, 1875 (2009).
  * Gilbert [1955] T. L. Gilbert, A lagrangian formulation of the gyromagnetic equation of the magnetization field, Phys. Rev. 100, 1243 (1955).
  * Gunnink  _et al._ [2022] P. M. Gunnink, B. Flebus, H. M. Hurst, and R. A. Duine, Nonlinear dynamics of the non-hermitian su-schrieffer-heeger model, Phys. Rev. B 105, 104433 (2022).
  * Hasan and Kane [2010] M. Z. Hasan and C. L. Kane, Colloquium: topological insulators, Reviews of modern physics 82, 3045 (2010).
  * Obana  _et al._ [2019] D. Obana, F. Liu, and K. Wakabayashi, Topological edge states in the su-schrieffer-heeger model, Phys. Rev. B 100, 075437 (2019).
  * Liu  _et al._ [2019b] S. Liu, W. Gao, Q. Zhang, S. Ma, L. Zhang, C. Liu, Y. J. Xiang, T. J. Cui, and S. Zhang, Topologically protected edge state in two-dimensional su–schrieffer–heeger circuit, Research 2019 (2019b).
  * Yuce and Ramezani [2019] C. Yuce and H. Ramezani, Topological states in a non-hermitian two-dimensional su-schrieffer-heeger model, Physical Review A 100, 032102 (2019).
  * Berger [1996] L. Berger, Emission of spin waves by a magnetic multilayer traversed by a current, Physical Review B 54, 9353 (1996).
  * Slonczewski [1996] J. C. Slonczewski, Current-driven excitation of magnetic multilayers, Journal of Magnetism and Magnetic Materials 159, L1 (1996).
  * Auerbach [2012] A. Auerbach,  _Interacting electrons and quantum magnetism_ (Springer Science & Business Media, 2012).
  * Abderraman Marrero and Aiat Hadj [2020] J. Abderraman Marrero and D. Aiat Hadj, Improving formulas for the eigenvalues of finite block-toeplitz tridiagonal matrices, Applied Mathematics and Computation 382, 125324 (2020).
  * Jin and Song [2019] L. Jin and Z. Song, Bulk-boundary correspondence in a non-hermitian system in one dimension with chiral inversion symmetry, Physical Review B 99, 081103 (2019).
  * Bender  _et al._ [1999] C. M. Bender, S. Boettcher, and P. N. Meisinger, $\mathcal{P}\mathcal{T}$-symmetric quantum mechanics, Journal of Mathematical Physics 40, 2201 (1999).