# Light Hybrid Meson Mixing and Phenomenology

E.S. Swanson  Department of Physics and Astronomy, University of Pittsburgh,
Pittsburgh, PA 15260, USA.

###### Abstract

A simple constituent model of gluodynamics that is motivated by lattice field
theory and the QCD Hamiltonian in Coulomb gauge is applied to descriptions of
hybrid meson flavor mixing and vector hybrid configuration mixing. Good
agreement with lattice gauge computations is obtained for flavor multiplet
masses, while mixing angles are in approximate agreement, given large errors.
The configuration mixing results are also in rough agreement with lattice
NRQCD calculations. Thus the viability of constituent gluon models of hybrid
hadrons and glueballs is supported. The results suggest that a flavor
multiplet of vector hybrids should appear with masses of approximately 2100,
2200, and 2300 MeV and that the isovector vector hybrid decay constant is
about 20 MeV. Similarly, the $\pi_{1}$ exotic hybrid should have isospin
partner states near 1760 and 1900 MeV, and it is suggested that the recently
seen $\eta_{1}$ hybrid signal is the latter state.

##  I Introduction

Although the notion of hadrons with explicit gluonic degrees of freedom has
been accepted for nearly 50 years, remarkably little is known about these
particles. Nevertheless, it is hoped that steady progress in lattice field
theory coupled with new effort at the GlueX, PANDA, and BESIII experiments
will finally shed light on these enigmatic states.

Models of gluonic degrees of freedom have generally assumed that they comprise
collective string-like excitations or some sort of quasiparticleMeyer:2015eta
. Recently, the quasiparticle picture has received support from lattice field
theory, where a measurement of the low lying charmonium hybrid spectrum
strongly suggests that an axial quasigluon with an effective mass something
less than 1000 MeV can explain the observed systematicsDudek:2011bn . This
intriguing observation has revived interest in constituent gluon models,
wherein early workHM  has evolved into more sophisticated modelling that
builds on QCDSzczepaniak:2001rg ; Swanson:1998kx ; Guo:2008yz ;
Szczepaniak:2003mr . This modelling starts with the QCD Hamiltonian in Coulomb
gauge and constructs gluonic quasiparticles with an Ansatz that builds gluonic
correlations in the vacuum. It is expected that the resulting field theory
admits reliable Fock space truncations which greatly enhances the ability to
model and compute hadronic properties.

The purpose of this work is to examine the viability of a simple constituent
gluon model of hybrid properties that is based on the considerations just
given. This will be done by computing flavor mixing of light hybrid mesons and
configuration mixing in vector mesons. The calculation is also of interest
because the flavor mixing mechanism is very different from that for canonical
mesons since the quark pair is in a color octet state. In particular, the
leading order mechanism annihilates and creates quark pairs via coupling to
the instantaneous Coulomb gauge interaction, while the next order mechanism
couples hybrid mesons to low lying glueballs. In this work, these glueballs
are described with the same degrees of freedom and dynamics as employed for
hybrids, thereby testing consistency of the model. Model validation is
possible because a comprehensive computation of the light meson spectrum in
lattice QCD has been made – a computation that includes isoscalar and
isovector low lying hybrid mesons and their mixing anglesDudek:2013yja .

A motivation for the constituent gluon model employed here will be given in
the next section. Section III applies the model to a computation of the light
meson spectrum, which fixes parameters for the light hybrid spectrum, which
permits investigation of flavor mixing in Section IV. Comparison to lattice
gauge results is made in Section V. A simple extension to vector hybrid-vector
($q\bar{q}$) meson mixing is presented in VI. A examination of the
implications of the results on the light hybrid and meson spectrum is made in
Section VII and we conclude in Section VIII.

##  II Constituent Gluon Model

###  II.1 Model Construction

A model dynamics capable of describing the interaction of constituent quarks
and axial quasigluons has been developedSzczepaniak:2001rg  and has been
applied to glueballsSzczepaniak:2003mr , the gluelump and hybrid
spectraGuo:2008yz ; Szczepaniak:2006nx ; Guo:2007sm , and to heavy hybrid
decaysFarina:2020slb . We briefly summarize the salient features of the model
to place the subsequent development in context.

The starting point for the dynamical model is the Hamiltonian of QCD in
Coulomb gauge. This gauge choice is expedient for model building because
Gauss’s law has been resolved, all degrees of freedom are physical, and an
explicit interaction potential that operates between quarks and gluons
emerges. This “Coulomb” interaction is written as

$V_{C}=\frac{1}{2}\int d^{3}x\,d^{3}y\,{\cal J}^{-1/2}\rho^{A}({\bf x}){\cal
J}^{1/2}\hat{K}_{AB}({\bf x},{\bf y};{\bf A}){\cal J}^{1/2}\rho^{B}({\bf
y}){\cal J}^{-1/2},$ (1)

where the Faddeev-Popov determinant is written as ${\cal J}\equiv{\rm
det}({\mbox{\boldmath$\nabla$}}\cdot{\bf D})$ and ${\bf D}$ is the adjoint
covariant derivative, ${\bf
D}^{AB}\equiv\delta^{AB}{\mbox{\boldmath$\nabla$}}-gf^{ABC}{\bf A}^{C}$. The
color charge density is given by

$\rho^{A}({\bf x})=f^{ABC}{\bf A}^{B}({\bf x})\cdot{\bf\Pi}^{C}({\bf
x})+\psi^{{\dagger}}({\bf x})T^{A}\psi({\bf x}).$ (2)

The kernel of the Coulomb interaction can be formally written asChrist:1980ku

$\hat{K}^{AB}({\bf x},{\bf y};{\bf A})\equiv\langle{\bf
x},A|\frac{g}{{\mbox{\boldmath$\nabla$}}\cdot{\bf
D}}(-{\mbox{\boldmath$\nabla$}}^{2})\frac{g}{{\mbox{\boldmath$\nabla$}}\cdot{\bf
D}}|{\bf y},B\rangle.$ (3)

Finally, ${\bf A}$ is the vector potential and ${\bm{\Pi}}$ is the conjugate
momentum given by the negative of the transverse chromoelectric field.

The Coulomb interaction, along with the quark and gluon kinetic energies,
gluon self-interactions, and the quark-transverse gluon interaction, $-g\int
d^{3}x\psi^{\dagger}\bm{\alpha}\cdot{\bf A}\psi$, comprise a full field-
theoretic version of QCD, with its accompanying nonperturbative features.

A quasigluon that is consistent with the constraints of QCD can be developed
by making a mean field model of the gluonic vacuum. The ensuing Schwinger-
Dyson equations can be truncated and solved to obtain estimates for the vacuum
expectation of the kernel, $\hat{K}^{AB}$, and for the gluon dispersion
relationshipSzczepaniak:2001rg . Here we choose to accept standard constituent
quark model phenomenology and lattice results for the static quark
interaction, and model the vacuum expectation of the Coulomb kernel as a
confining potential:

$\langle\hat{K}^{AB}({\bf r},{\bf
A})\rangle\to\delta^{AB}\left(-\frac{3}{4}\mathcal{C}+\frac{a_{S}}{r}-\frac{3}{4}\sigma
r\right).$ (4)

Of course this reproduces the successes of the Cornell potential in
nonrelativistic quark models. Higher terms in the $n$-body expansion of
$\hat{K}$ can be incorporated in the formalism as required.

The vacuum model also gives rise to a quasigluon that can be described by a
field expansion parameterized with a dispersion relationship,
$\omega=\omega(k)$. Direct computation in the vacuum Ansatz yields an
expression that is well-approximated bySzczepaniak:2001rg

$\omega^{2}=k^{2}+m_{g}^{2}{\rm e}^{-k/b_{g}}$ (5)

where the dynamical gluon mass is $m_{g}\approx 600$ MeV and the parameter
$b_{g}\approx 6000$ MeV. We stress that the gluon remains transverse and
properties, such as Yang’s theorem, remain in place. Other vacuum Ansätze are
possible, for example a Gaussian wavefunctional (equivalent to the mean field
approximation described) can be combined with the Faddeev-Popov operator,
which gives rise to a dispersion relation that is well-described by the Gribov
form, $\omega^{2}=k^{2}+m_{g}^{4}/k^{2}$Feuchter:2004mk .

The model can be validated by computing the excited gluonic potentials in the
case of fixed color sources. Doing so reveals that the potential surfaces are
not ordered according to lattice results. The discrepancy can be corrected by
including tri-linear gluonic terms in $\hat{K}^{AB}$ in the
computationSzczepaniak:2006nx . These contributions are zero for the lowest
surface, which are dealt with exclusively in this work, thus trilinear
couplings are neglected.

The resulting model can be thought of as a minimal extension of the
constituent quark model with the addition of constituent gluon degrees of
freedom and possible additional couplings (for example, the tri-linear gluon
coupling or the gluon–Coulomb interaction).

###  II.2 Hybrid States

As is traditional, it is assumed that hybrid mesons are dominated by Fock
states with the lowest number of constituents. This approximation is
unreasonable in perturbative QCD but is made plausible here by the relatively
large quasigluon mass. As stated above, one of the goals of this study is to
test this statement.

It is convenient to construct the total gluon spin, $j_{g}$, by coupling the
gluon spin projection to the gluon angular momentum, $\ell_{g}$. Converting to
the gluon helicity basis and assuming that $\ell_{g}=j_{g}$ reduces the
product of two Wigner matrices to one and produces a factor of

$\chi^{(-)}_{\lambda,\mu}\equiv\langle
1\lambda\ell_{g}0|\ell_{g}\mu\rangle=\begin{cases}0,\ell_{g}=0\\\
\frac{\lambda}{\sqrt{2}}\delta_{\lambda,\mu},\ell_{g}\geq 1\end{cases}.$ (6)

This represents a transverse electric (TE) gluon in our model and forms the
explicit realization of the axial constituent gluon. Alternatively, one may
set $\ell_{g}=j_{g}\pm 1$ and obtain a transverse magnetic (TM) gluon with a
Clebsch factor given by
$\chi^{(+)}_{\lambda,\mu}=\delta_{\lambda,\mu}/\sqrt{2}$. Here we work
exclusively with low lying TE hybrid mesons.

Combining with quark spins yields the final expression for a hybrid creation
operator

$\displaystyle|JM[LS\ell
j_{g}\xi]\rangle=\frac{1}{2}T_{ij}^{A}\int\frac{d^{3}q}{(2\pi)^{3}}\,\frac{d^{3}k}{(2\pi)^{3}}\,\Psi_{j_{g};\ell
m_{\ell}}({\bf k},{\bf
q})\,\sqrt{\frac{2j_{g}+1}{4\pi}}\,D_{m_{g}\mu}^{j_{g}*}(\hat{k})\,\chi^{(\xi)}_{\mu,\lambda}$
$\displaystyle\times\langle\frac{1}{2}m\frac{1}{2}\bar{m}|SM_{S}\rangle\,\langle\ell
m_{\ell},j_{g}m_{g}|LM_{L}\rangle\,\langle SM_{S},LM_{L}|JM\rangle\,b_{{\bf
q}-\frac{{\bf k}}{2},i,m}^{\dagger}\,d_{-{\bf q}-\frac{{\bf
k}}{2},j,\bar{m}}^{\dagger}\,a^{\dagger}_{{\bf k},A,\lambda}|0\rangle.$ (7)

Quark and gluon particle operators are understood to create quasiparticles and
$|0\rangle$ refers to the correlated vacuum discussed above.

By construction, the hybrid state is an eigenstate of parity and charge
conjugation with eigenvalues given by

$P=\xi(-1)^{\ell+j_{g}+1}\ \ \textrm{and}\ \ C=(-1)^{\ell+S+1}.$ (8)

###  II.3 Glueball States

A reasonably large quasigluon mass encourages modelling charge-conjugation
positive (negative) glueballs as two (three) quasigluon states. Again, it is
preferable to work in the helicity basis, where much of the algebra
simplifies. Combining two gluons into states with good parity and total
angular momentum $J$ can be achieved withSzczepaniak:2003mr

$|JM;\eta\rangle=\frac{1}{\sqrt{2}}\left(|JM;\lambda,\lambda^{\prime}\rangle+\eta|JM;-\lambda,-\lambda^{\prime}\rangle\right),$
(9)

for which $P|JM;\eta\rangle=\eta(-)^{J}\,|JM;\eta\rangle$ and $\eta=\pm 1$.

The helicity states are constructed as

$|JM;\lambda,\lambda^{\prime}\rangle=\frac{1}{\sqrt{2(N_{c}^{2}-1)}}\sqrt{\frac{2J+1}{4\pi}}\,\int\frac{d^{3}k}{(2\pi)^{3}}\,\psi(k)\,D^{J*}_{M,\lambda-\lambda^{\prime}}(\phi,\theta,-\phi)\,\Pi\,a^{\dagger}({\bf
k},\lambda,A)a^{\dagger}(-{\bf k},\lambda,A)|0\rangle.$ (10)

The number of colors is denoted $N_{c}$, $A$ is an adjoint color index, and
$\Pi$ is a Jacob-Wick phase that will not be important to the following
development. The wavefunction $\psi$ is determined by solving the Tamm-Dancoff
equation that is obtained by evaluating the QCD Hamiltonian in the appropriate
$J^{P}$ channel. This yields the leading order contribution involving the
Coulomb interaction, Eq. 1. Higher order contributions from gluon exchange and
the four-gluon interaction can be incorporated if desired. The resulting
spectrum is reported in Ref. Szczepaniak:2003mr , where it is compared to
lattice field theory computations. The spectra agree quite well where they
overlap, with the largest deviation being about 200 MeV.

##  III Light Hybrid Spectrum

###  III.1 Model Parameter Selection

Because the primary goal is to compute hybrid mixing masses and angles, it is
not necessary to obtain a precise hybrid spectrum. This is convenient because
very little is known experimentally and because lattice field computations at
physical pion masses and with coupled channel effects are not yet available.
We therefore focus on spin-averaged hybrid multiplets in the following. In
practice this means neglecting transverse gluon exchange contributions to the
quark and gluon interactions. Thus the model parameters are the quark mass,
the Coulomb coefficient, the string tension, and the constant shift: $m$,
$a_{S}$, $\sigma$, $\mathcal{C}$. Recall that the gluonic parameters $m_{g}$
and $b_{g}$ have been fixed by the vacuum model. Including transverse gluons
will introduce the coupling $g$ as well, which can be fixed by
$a_{S}=g^{2}/(4\pi)$ or from other considerations to be discussed.

Model parameters will be fixed by fitting to 16 light isovector nonexotic
meson masses. As a check of stability we also fit to 30 isospin 0, 1/2, and 1
light mesons whose identities are reasonably well established. This is not
necessarily a simple procedure since identifying “non-canonical” properties in
the light mesons is notoriously difficult. Famous examples include the
$f_{0}(500)$ which has come and gone in the Review of Particle Physics (PDG)
over the years. Similarly, the $a_{0}(980)$ has been identified as a
$q\bar{q}$ state, a tetraquark, or a $K\bar{K}$ bound state by various
authors. Lastly, the pion is a pseudo-Goldstone boson, hence simple quark
models cannot be expected to reproduce its properties.

We do not presume to have a definitive description of light mesons and
therefore will fit several model variations to obtain a sense of parameter
stability in the subsequent work. These models are (i) spin-independent
interaction (Eq. 4) with a smeared hyperfine interaction, (ii) spin-
independent interaction, (iii) spin-dependent interaction, (iv) variation on
(iii), (v) spin-dependent interaction fit to 30 mesons. For the sake of
comparison, results for models with relativistic quark kinetic energies are
also given below, although these are not used in the subsequent analysis.

The hyperfine interaction used in model (i) is given by

$V_{H}=\frac{32\pi
a_{S}}{9m_{q}m_{\bar{q}}}\left(\frac{\varsigma}{\sqrt{\pi}}\right)^{3}\,\exp(-\varsigma^{2}r^{2})\,{\bf
S}_{q}\cdot{\bf S}_{\bar{q}},$ (11)

with the smearing parameter set to $\varsigma=0.897$ GeV. For the other
models, the spin-dependent interaction is defined by
$V_{SD}=V_{H}+V_{LS}+V_{T}$ with:

$V_{H}=\frac{8\alpha_{H}}{3}b_{0}^{2}\frac{\textrm{e}^{-b_{0}r}}{3m_{q}m_{\bar{q}}r}{\bf
S}_{q}\cdot{\bf S}_{\bar{q}},$ (12)

$V_{LS}=\left(\frac{4\alpha_{H}}{3\rho^{3}}+\frac{\epsilon\sigma}{\rho}\right)\frac{{\bf
L}\cdot{\bf
S}}{m_{q}m_{\bar{q}}}+\frac{1}{2}\left(\frac{4}{3}\frac{\alpha_{H}}{\rho^{3}}+\frac{(2\epsilon-1)\sigma}{\rho}\right)\,\left(\frac{{\bf
L}\cdot{\bf S}_{q}}{m_{q}^{2}}+\frac{{\bf L}\cdot{\bf
S}_{\bar{q}}}{m_{\bar{q}}^{2}}\right),$ (13)

and

$V_{T}=\frac{4\alpha_{H}}{3m_{q}m_{\bar{q}}\rho^{3}}\left({\bf
S}_{q}\cdot\hat{r}\,{\bf S}_{\bar{q}}\cdot\hat{r}-{\bf S}_{q}\cdot{\bf
S}_{\bar{q}}\right).$ (14)

The ultraviolet singularity in these expressions has been regulated by
freezing $r$ at $r_{0}$ once $r<r_{0}$; this is denoted as $\rho$ in the
equations. The parameter $\epsilon$ that appears in the spin-orbit tensor
interaction represents a mixture of “scalar” and “vector” confinement models.
Lattice computations find that $\epsilon\approx 0.25$KK , which is used in
model (v).

Results for the fits are presented in Table 1 and will be used to model hybrid
mesons.

Table 1: Model Parameters for the Isovector Meson Spectrum. model |  $m\ (m_{s})$ (MeV) | $a_{S}$ |  $\sigma$ (GeV2) |  $\mathcal{C}$ (MeV) | $\alpha_{H}$ |  $b_{0}$ (GeV-1) |  $r_{0}$ (GeV-1) |  $\epsilon$ (GeV-1) | rel error | avg deviation (MeV)
---|---|---|---|---|---|---|---|---|---|---
i. [hyp] | 335 | 0.59 | 0.16 | -697 | – | – | – | 0 | 9% | 94
ii. [SI] | 300 | 1.52 | 0.071 | 110 | – | – | – | 0 | 7% | 66
iii. [SD] | 400 | 1.8 | 0.06 | 230 | 1.3 | 0.60 | 4.7 | 0 | 5% | 54
iv. [SD] | 330 | 2.1 | 0.055 | 385 | 1.3 | 0.53 | 5.2 | 0 | 5% | 54
v. [30] | 375 (525) | 1.3 | 0.059 | 112 | 0.77 | 0.84 | 5.2 | 0.25 | 7% | 77
Rel/SI | 200 | 0.59 | 0.14 | -246 | – | – | – | 0 | 6% | 59
Rel/SD | 400 | 0.72 | 0.14 | -359 | 1.1 | 0.20 | 4.4 | 0 | 5% | 55

###  III.2 Hybrid Mesons

Spin-independent hybrid wavefunctions are obtained by considering the
nonrelativistic limit of the interaction of Eq. 1 (with Eq. 4). The resulting
spectrum can be categorized according to interpolating operators, as indicated
in Table 2JKM . Here ${\bf B}$ is the chromomagnetic field, and $\psi$ and
$\chi$ are heavy quark and antiquark fields, respectively. The remaining
columns give the corresponding quantum numbers in the present model and the
hybrid meson quantum numbers in the specified multiplet.

Table 2: $J^{PC}$ Hybrid Multiplets. multiplet | operator | $\xi$ | $j_{g}$ | $\ell$ | $L$ | $J^{PC}\ S=0\ (S=1)$
---|---|---|---|---|---|---
$H_{1}$ | $\psi^{\dagger}{\bf B}\chi$ | -1 | 1 | 0 | 1 |  $1^{--}$, $(0,1,2)^{-+}$
$H_{2}$ | $\psi^{\dagger}{\mbox{\boldmath$\nabla$}}\times{\bf B}\chi$ | -1 | 1 | 1 | 1 |  $1^{++}$, $(0,1,2)^{+-}$
$H_{3}$ | $\psi^{\dagger}{\mbox{\boldmath$\nabla$}}\cdot{\bf B}\chi$ | -1 | 1 | 1 | 0 |  $0^{++}$, $(1^{+-})$
$H_{4}$ | $\psi^{\dagger}[{\mbox{\boldmath$\nabla$}}{\bf B}]_{2}\chi$ | -1 | 1 | 1 | 2 |  $2^{++}$, $(1,2,3)^{+-}$

The bound state equation is obtained from the model QCD Hamiltonian by
computing the expectation value, $\langle
J^{\prime}M^{\prime}[L^{\prime}S^{\prime}\ell^{\prime}j^{\prime}\xi^{\prime}]|H|JM[LS\ell
j\xi]\rangle$. As mentioned, we seek spin-independent multiplets and therefore
consider the nonrelativistic limit of the currents in Eq. 1. This gives rise
to instantaneous quark-antiquark and (anti)quark-gluon interactions that
generate the bound state.

A novel method for solving the quantum mechanical three-body problem was
applied to solve the resulting Schrödinger equation. This consisted of writing
the hybrid wavefunction as a sum over a product Ansatz of the form

$\Psi_{j_{g};\ell m_{\ell}}({\bf k},{\bf
q})=\chi_{j_{g}}(k)\,\varphi_{\ell}(q)Y_{\ell,m_{\ell}}(\hat{q}).$ (15)

This makes explicit the angular momentum dependence in the $q$ coordinate,
while the gluon angular momentum dependence is contained in the Wigner
rotation matrix in Eq. 7. In practice the basis used is nearly diagonal in the
quantum numbers, having only a coupling between the $H_{1}$ and $H_{2}$
multiplets induced by mixing between TM and TE hybrids (for spin-independent
interactions). This will be small and is neglected. It is thus possible to
solve for $\varphi$ and $\chi$ separately and iterate the coupled equations to
convergence.

We sketch the idea here, ignoring all indices for simplicity. Write the
Hamiltonian generically as $K_{q}+K_{g}+V$, where the first two terms are the
quark and gluon kinetic energy operators obtained from the QCD Hamiltonian and
the potential includes the sum over the three possible instantaneous
interactions. Then vary
$\langle\varphi\chi|H|\varphi\chi\rangle+\lambda(\langle\varphi\chi|\varphi\chi\rangle-1)$
with respect to $\varphi$ and $\chi$. Eliminating the Lagrange multiplier
yields

$\displaystyle
K_{q}\varphi+\int\chi^{*}K_{g}\chi\cdot\varphi+\int\chi^{*}V\chi\cdot\varphi$
$\displaystyle=E\varphi$ $\displaystyle
K_{g}\chi+\int\varphi^{*}K_{q}\varphi\cdot\chi+\int\varphi^{*}V\varphi\cdot\chi$
$\displaystyle=E\chi.$ (16)

Eq. 16 reduces to two one-dimensional equations. We solve this system by using
a simple and accurate discretization methodBS , diagonalizing the Laplacian
operator to deal with the (momentum space) gluonic kinetic energyBS , and
iterating. The latter step requires an initial guess for $\varphi$ and $\chi$,
which is obtained variationally. In practice the method converges very
quickly, and a precise solution to the quantum bound state problem is
obtained.

Calculations were done with the gluon dispersion relationship of Eq. 5 with
$m_{g}$ set to 600 MeV. Performing the same calculations with the Gribov form
of the dispersion relationship and the same mass scale yielded very similar
hybrid masses, with typical results being approximately 10 MeV heavier than
given with Eq. 5.

The resulting hybrid multiplet mass splittings with respect to $H_{1}$ are
shown in Fig. 1. Detailed model validation is not feasible at present. A
lattice computation of the light meson spectrum at a pion mass of 391 MeV
exists, but it has not been able to distinguish enough hybrid states to
determine the spin-averaged spectrumDudek:2013yja . This calculation does,
however, set the scale for the spin-averaged $H_{1}$ multiplet (taken to be
the $S=0$ $1^{--}$ mass) to be $2190\pm 20$ MeV.

On the experimental side, the situation is even more sparse and confused. Past
claims to exotic $\pi_{1}$ states near 1400 and 1600 MeVWorkman:2022ynf  have
recently been challenged, with a consensus emerging that only one $\pi_{1}$
exists near 1600 MeV (results from two analyses are $M=1564\pm 24\pm 86$ MeV,
$\Gamma=492\pm 54\pm 102$ MeV JPAC:2018zyd  and $M=1623\pm 47\pm 50$ MeV,
$\Gamma=455\pm 88\pm 150$ MeV Kopf:2020yoa ).

Figure 1 displays a very large predicted splittings, $H_{2}-H_{1}\approx 500$
MeV. This splitting can be estimated from the lattice calculation of Ref.
Dudek:2013yja  using the $J^{PC}=2^{+-}$ to $1^{--}$ mass difference (assuming
that the $2^{+-}$ mass is approximately the $1^{+-}$ mass, which is supported
by lattice calculations at the charmonium massHadronSpectrum:2012gic ;
Cheung:2016bym ), yielding a value of approximately 250 MeV. Notice also that
$H_{3}$ lies above $H_{4}$. This situation also occurs in a similar model
calculation for charmonium hybrids, reported in Ref. Farina:2020slb , where it
is seen to disagree with the lattice ordering, $H_{4}>H_{3}$. Thus, although
the model broadly agrees with lattice field theory calculations, where
available, it is appears that additional effects, such as occur at higher
order in the $1/m$ expansion, and further model tuning may be important to
obtaining detailed agreement.

Figure 1: Hybrid Multiplet Mass Splittings and $H_{1}$ Variational Parameters.
Model i (red), ii (black), iii (green), iv (blue),v (purple).

It is useful to find Gaussian estimates to the exact wavefunctions considered
here so that an analytic evaluation of the mixing matrix elements (to be
discussed in the next section) can be made. For this purpose we write
$\chi_{j_{g}=1}(k)\propto k\exp(-k^{2}/\beta_{g}^{2})$ and
$\varphi_{\ell}(q)\propto q^{\ell}\exp(-q^{2}/\alpha_{q}^{2})$. The parameters
are estimated by optimizing the energy variationally and are shown to the
right in Fig. 1 for the five models considered here.

##  IV Hybrid Flavor Mixing

The topic of light meson flavor mixing is replete with experimental data, with
much to be gleaned from a plethora of decay modes. In contrast, theoretical
understanding of the issue is essentially nonexistent. The only certainty is
that mixing occurs via nonperturbative gluodynamics, presumably dominated by
coupling to intermediate glueballs. Since so little is known of glueballs and
their dynamics, theory is largely guesswork. (We remark that mixing in the
$\eta-\eta^{\prime}$ system is unique in that the axial anomaly makes a
quantifiable nonperturbative contribution.)

In contrast to the rather grim situation with canonical mesons, hybrid mesons
have their dominant quark configuration in a color octet, and can therefore
mix perturbatively. The general situation in a given isospin multiplet
involves amplitudes for mixing $u\bar{u}\leftrightarrow u\bar{u}$ and
$u\bar{u}\leftrightarrow d\bar{d}$. These are expected to be nearly identical,
hence both are labelled $A_{nn}$. Mixing $u\bar{u}\leftrightarrow s\bar{s}$
will be denoted $A_{ns}$, while $s\bar{s}\leftrightarrow s\bar{s}$ will be
$A_{ss}$.

As we have mentioned, mixing to positive charge conjugation glueballs is first
order in the strong coupling, and can therefore expected to be important. We
label these amplitudes $\mathcal{A}^{(n)}_{f}$ where $f=n,s$ denotes the
annihilated quark flavor and $(n)$ denotes a radially excited glueball of the
relevant quantum numbers. We will show that this sum saturates quickly, so
only the ground state glueball is considered in the following.

Diagonal elements of the QCD Hamiltonian will be written as $m$ for the
$u\bar{u}$ and $d\bar{d}$ cases, $m+\Delta m$ for $s\bar{s}$, and $M_{gb}$ for
the bare glueball mass. Thus, in the $u\bar{u}$, $d\bar{d}$, $s\bar{s}$ basis
the matrix elements of the QCD Hamiltonian are

$H_{uds}=\begin{pmatrix}m+A_{nn}&A_{nn}&A_{ns}&\mathcal{A}^{(0)}_{n}\\\
A_{nn}&m+A_{nn}&A_{ns}&\mathcal{A}^{(0)}_{n}\\\ A_{ns}&A_{ss}&m+\Delta
m+A_{ss}&\mathcal{A}^{(0)}_{s}\\\
\mathcal{A}^{(0)}_{n}&\mathcal{A}^{(0)}_{n}&\mathcal{A}^{(0)}_{s}&M_{gb}\end{pmatrix}.$
(17)

Switching to the isospin basis $(u\bar{u}-d\bar{d})/\sqrt{2}$,
$(u\bar{u}+d\bar{d})/\sqrt{2}$, partially diagonalizes the mass matrix:

$H_{iso}=\begin{pmatrix}m&0&0&0\\\
0&m+2A_{nn}&\sqrt{2}A_{ns}&\sqrt{2}\mathcal{A}^{(0)}_{n}\\\
0&\sqrt{2}A_{ns}&m+\Delta m+A_{ss}&\mathcal{A}^{(0)}_{s}\\\
0&\sqrt{2}\mathcal{A}^{(0)}_{n}&\mathcal{A}^{(0)}_{s}&M_{gb}\end{pmatrix}.$
(18)

A final diagonalization then gives the isovector, isoscalar, and glueball
masses and mixing angles.

The leading (order $g^{2}$) sources of hybrid mixing are direct annihilation
(Fig. 2, left), which occurs at first order in perturbation theory, second
order mixing via glueball states (Fig. 2, right), or second order mixing to
the meson-meson continuum. Flavor mixing via coupling to the meson-meson
continuum has been enigmatic since the beginnings of the quark model. The
issue, first stressed by LipkinZvi , is that continuum mixing can vitiate the
Okubo-Zweig-Iizuka (OZI) rule because a process such as $J/\psi\to
D\bar{D}\to\omega$ is not suppressed. How the OZI rule arises in spite of this
mechanism has been explored by Geiger and Isgur, who argue that cancellations
occur when all possible intermediate meson-meson channels are summed, giving
rise to an emergent scale that is much smaller than
$\Lambda_{QCD}$Geiger:1991ab ; Geiger:1992va . We shall assume that hybrid
mixing via coupling to the meson-meson continuum is similarly suppressed in
the case of hybrid states.

The leading order expression for the mixing amplitude, shown in Fig. 2, is
given by

$A_{ff^{\prime}}=\frac{1}{m_{f}m_{f^{\prime}}}\int\frac{k^{2}dk}{(2\pi)^{3}}\,\frac{d^{3}q}{(2\pi)^{3}}\,\frac{d^{3}q^{\prime}}{(2\pi)^{3}}\,\Psi_{f}(k,{\bf
q})\Psi_{f^{\prime}}^{*}(k,{{\bf
q}}^{\prime})k^{2}V(k)\,B_{J}=\frac{F_{f}F_{f^{\prime}}}{8m_{f}m_{f^{\prime}}}\int\frac{k^{2}dk}{(2\pi)^{3}}\,|\chi_{1}(k)|^{2}k^{2}V(k)\,B_{J}.$
(19)

The potential $V(k)$ is the Fourier transform of Eq. 4. Wigner rotation
matrices have been integrated and Clebsch-Gordan sums have been done to give
the first form. The second follows from the product Ansatz of Eq. 15 and
introduces the “octet decay constant”
$F_{f}=\int\frac{d^{3}q}{(2\pi)^{3}}\varphi_{\ell=0}({\bf q})$, where implicit
flavor-dependence is labelled by $f$. Evaluation of the discrete sums is
considerably simplified because the quark vertex forces $S=S^{\prime}=1$ and
$\ell=\ell^{\prime}=0$. Thus hybrid mixing at this order only exists in the
spin-triplet portion of the $H_{1}$ multiplet with relative strengths given by
the Clebsch factor

$B_{J}=\begin{cases}0&J=0\\\ 1&J=1\\\ 3/5&J=2\end{cases}.$ (20)

Finally, the integrals can be performed if Gaussian approximate wavefunctions
are employed, giving

$A_{ff^{\prime}}=\frac{\pi
B_{J}}{2m_{f}m_{f^{\prime}}}\left(\frac{\alpha_{f}\alpha_{f^{\prime}}}{\pi}\right)^{3/2}[a_{S}+\frac{b}{\beta_{g}^{2}}],$
(21)

where $\beta_{g}$ is the gluonic scale introduced after Eq. 16.

Figure 2: (left) First order mixing diagram. (right) Hybrid-gluball mixing
(crossed diagram not shown).

Second order mixing via intermediate glueballs can be computed with the
amplitude of Fig. 2 (right). This diagram also features quark-antiquark
annihilation in the ${}^{3}S_{1}$ channel, and therefore mixing is predicted
to be zero for hybrids in the $H_{2}$, $H_{3}$, and $H_{4}$ multiplets, as
well as the light vector, $H_{1}(1^{--})$ multiplet. This striking observation
can only be reasonably evaded by mixing to negative charge conjugation (three-
quasigluon) glueballs, which is expected to be small, or by coupling to the
meson-meson continuum, which remains an enigmatic feature of hadronic physics.

The expression for the amplitude coupling a hybrid to the $n$th radial
glueball excitation is

$\mathcal{A}^{(n)}_{f}=-\frac{i[gF_{f}]}{4}\,\int\frac{k^{2}dk}{(2\pi)^{3}}\,\frac{\psi_{n}^{*}(k)\chi(k)}{\sqrt{\omega(k)}}\,C_{J},$
(22)

where

$C_{J}=\begin{cases}4,&J=0\\\ 0,&J=1\\\ 4/\sqrt{10},&J=2\end{cases}.$ (23)

Note that the octet decay constant, $F_{f}$, appears in the glueball
amplitude, this time combined with the strong coupling constant. Evidently,
obtaining an accurate estimate of the decay constant is important. This can be
problematic because it is known that the nonrelativistic approximation over-
estimates the value of (traditional) meson decay constants. Incorporating
relativistic effects helps, but it appears that further softening is required.
This softening can occur, for example, via the effect of the running coupling
on the wavefunction at the originLakhina:2006vg . These issues are exacerbated
in the case of glueball mixing (Eq. 22) because the strong coupling constant
is explicit. One can set the value of the strong coupling from the model via
$g=\sqrt{4\pi a_{S}}$. This has the appeal of consistency, but might not be
optimal because $a_{S}$ is a model parameter that is determined by bulk light
hadron properties that are dominated by distances near 1 fm. Alternatively, a
decay constant is a short range phenomenon–roughly speaking we wish to
evaluate $[gF_{f}]\sim\sqrt{\alpha_{S}(r=0)}\varphi(r=0)$. Of course this
implies a nontrivial infrared fixed point for the running coupling, which will
be assumed here. A better way to proceed is to write
$[gF_{f}]\sim\alpha_{S}(Q_{*})\varphi(r=0)$ where $Q_{*}$ is a scale that is
tuned to the physics. An alternative, that is adopted here, is to implicitly
fix the scale by averaging; thus we set

$[gF_{f}]=\int\frac{d^{3}q}{(2\pi)^{3}}\,\sqrt{4\pi\alpha_{V}(q)}\,\phi_{\ell=0}(q)$
(24)

where the running coupling is parameterized as

$\alpha_{V}(q)=\frac{4\pi}{b_{0}\log(\frac{q^{2}+M^{2}}{\Lambda_{V}^{2}})}.$
(25)

with $b_{0}=11-2n_{f}/3=9$ for our case. This is a reasonably common model
that has been advocated for renormalizing exclusive processesBrodsky:1997dh
and has been used in modelling heavy meson propertiesLakhina:2006vg .
Parameters chosen in these studies were $M=870$ MeV, $\Lambda_{V}=160$ MeV or
$M=1000$ MeV, $\Lambda_{V}=250$ MeV, respectively. Results for both model
choices ($\sqrt{4\pi a_{S}}$ and $\sqrt{4\pi\alpha_{V}}$) will be presented in
the following section.

Lastly, we address the issue of the convergence of the hybrid mixing amplitude
in the sum over glueball excitations. The perturbative glueball mixing
amplitude is given by

$A^{gb}_{ff^{\prime}}=\sum_{n}\frac{{\mathcal{A}_{f}^{(n)}}^{*}\mathcal{A}^{(n)}_{f^{\prime}}}{M_{hyb}-M^{(n)}_{gb}}.$
(26)

This sum is expected to converge quickly because the integral in Eq. 22
rapidly decreases with radial quantum number. In fact, the integral would be
zero (for $n>0$) in a simple harmonic approximation to the wavefunctions if
the glueball and hybrid scales were the same. An explicit calculation in the
$J^{PC}=0^{-+}$ case shows that $\mathcal{A}_{n}^{(1)}\approx
0.37\mathcal{A}_{n}^{(0)}$. The corresponding term in $A^{gb}_{nn}$ is further
suppressed by the larger radial glueball mass, giving a final contribution
that is only 1% of the leading term in the sum. Because of this we only
considered hybrid coupling to the lowest mass glueball in a given $J^{P}$
channel in Eq. 18.

##  V Comparison to Lattice Computations

We proceed by diagonalizing the matrix of Eq. 18 for the three nontrivial
cases, $J^{PC}=0^{-+}$, $1^{-+}$, and $2^{-+}$ (recall that all other hybrid
mesons are predicted to have negligible mixing). The resulting masses are
shown in Fig. 3, along with lattice masses computed in Ref. Dudek:2013yja .
The latter are computed at a pion mass of 391 MeV and therefore may experience
some shifts in going to physical quark masses. In view of this, the model
results have been normalized by setting $m$ to the isovector mass (shown in
blue) for each $J^{PC}$ multiplet. Model isoscalar masses are also shown as
black oblongs in the figure (states dominated by glueball components lie
substantially higher and are not shown). Computing masses over the models of
Table 1 gives an indication of parameter dependence. This dependence is
indicated in the figure by vertical grey bars.

Figure 3: Isovector and isoscalar hybrid masses. Lattice results are blue
(isovector) and gray boxes (isoscalars). The box heights indicate statistical
uncertainty. Model results for isoscalar masses are shown in black with model
variation as vertical bars. Stars indicate model results with the running
octet decay constant model of Eq. 24.

As discussed above, no mixing is expected outside the $H_{1}$ multiplet.
Evidently this prediction agrees very well with the lattice results for the
$J^{PC}=0^{+-}$, $2^{+-}$, and $2^{+-^{\prime}}$ multiplets. Countering this
is the $1^{--}$ multiplet, where the lowest isoscalar is computed to be
approximately 140 MeV above the isoscalar, rather than degenerate with it as
predicted here. This curious situation is difficult to reconcile with the
current model. The dominant mixing effect would be via negative parity
glueballs, which requires a gluon emission followed by a quark spin flip and
then quark pair annihilation. This process will be suppressed by a relatively
large three-quasigluon glueball mass and the spin flip.

Turning attention to the $2^{-+}$ multiplet, reasonable agreement is seen for
the lowest isoscalar mass, with some overlap for the higher
isoscalar–especially for the preferred running coupling octet decay constant
of Eq. 24. In this case we suspect that the higher lattice mass is anomalously
low since it is implausible for isoscalars to shift mass with respect to the
isovector while remaining degenerate. In contrast, the $1^{-+}$ multiplet
splits as expected and in reasonable agreement with the model calculation.

Finally, the $0^{-+}$ multiplet is very unusual–the lattice results imply very
small mixing, in contrast the $\sqrt{a_{S}}$ model predicts that the light
isoscalar lies from 50 to several hundred MeV below the isovector. This is
because there is no direct mixing, $A_{nn}(J^{PC}=0^{-+})=0$, and mixing with
the pseudoscalar glueball drives the light isoscalar mass down. The upward
shift of the glueball is comparable to the downward shift of the light
isoscalar, hence if the larger isoscalar shift proves correct, then unquenched
lattice calculations of the pseudoscalar glueball mass should find it shifted
several hundred MeV above the unquenched value of approximately 2600
MeVChen:2005mg . Of course the lattice results of Fig. 3 argue against this,
and imply that the splittings are at the small end of the predicted range.
Alternatively, the preferred running octet decay constant model gives a much
smaller mixing that is reasonably close to lattice results–although still
indicating a novel light isoscalar.

$J^{PC}$ | nominal state | $u\bar{u}g$ | $s\bar{s}g$ | $gg$
---|---|---|---|---
$1^{--}$ | light | $\approx 100$ | $\approx 0$ | $\approx 0$
| heavy | $\approx 0$ | $\approx 100$ | $\approx 0$
| glueball | $\approx 0$ | $\approx 0$ | $\approx 100$
$0^{-+}$ | light | 62 [87] | 17 [6] | 21 [7]
| heavy | 24 [8] | 75 [91] | 0.5 [1]
| glueball | 13 [5] | 8 [3] | 78 [92]
$1^{-+}$ | light | 37 | 63 | 0
| heavy | 63 | 37 | 0
| glueball | 0 | 0 | 100
$2^{-+}$ | light | 54 [59] | 46 [41] | 0.1 [0]
| heavy | 32 [38] | 67 [61] | 1 [1]
| glueball | 2 [0.4] | 1 [0.6] | 97 [99]
Table 3: Model Hybrid $H_{1}$ Multiplet Fock Space Components (%). Results for
the running octet decay constant model of Eq. 24 are shown in square brackets.

The authors of Ref. Dudek:2013yja  also report mixing angles between the two
lightest isoscalar mesons. These were computed under the assumption that the
states do not mix with nearby glueballs. Of course this assumption is not made
here; however, the erstwhile mixing angle can still be computed. This will
lead to an ambiguous result if mixing to glueballs is substantial. As shown in
Tab. 3, we find that this is not the case for all multiplets except
$J^{PC}=0^{-+}$ (in the $\sqrt{a_{S}}$ model). Lattice mixing angles for three
volumes are shown in blue in Fig. 4; model results are displayed as black
oblongs. As with the meson masses, the results are in broad agreement with
lattice (in view of the large errors), with the largest discrepancy being in
the vector multiplet again.

Figure 4: Lattice and model mixing angles. Lattice mixing angles for three
volumes($16^{3}$, $20^{3}$, $24^{3}$) with statistical errors indicated by the
box height (blue). Model results (black) with model variation shown as
vertical bars. Stars indicate model results with the running octet decay
constant model of Eq. 24.

##  VI Hybrid-Canonical Meson Mixing

The encouraging results of the previous section motivate the consideration of
hybrid meson mixing with canonical mesons. Work of this sort dates back to the
beginnings of the bag model and QCD, see for exampleBarnes:1977hg . For a
computation based the Born-Oppenheimer approach, see Ref. Gerasimov:1998tm .

Following the philosophy advocated here, this process will be mediated by
gluon production from quark and anti-quark lines. Thus we seek
$\mathcal{H}\equiv\langle q\bar{q}g|ig\int\psi^{\dagger}\bm{\alpha}\cdot{\bf
A}\psi|q{\bar{q}}\rangle$. Configuration mixing of this sort is of most
interest for vector states since it has implications on the coupling of vector
hybrids to $e^{+}e^{-}$.

Taking the nonrelativistic limit of the vertex, performing integrals over the
angles of ${\bf k}$, and doing the Clebsch sums gives a result involving an
integral of the hybrid wavefunction convoluted with the vector quarkonium
wavefunction:

$\mathcal{H}=-i\frac{g}{m}\frac{2\sqrt{4\pi}}{3}\int\frac{d^{3}q}{(2\pi)^{3}}\,\frac{k^{2}dk}{(2\pi)^{3}}\,\frac{k}{\sqrt{\omega(k)}}\Psi^{*}(k,{\bf
q})\psi({\bf q}+{\bf k}/2).$ (27)

The hybrid wavefunction is obtained with the method of Sect. III.2 for
$u\bar{u}g$, $c\bar{c}g$, and $b\bar{b}g$ $H_{1}$ vector hybrids.
Wavefunctions corresponding to $\rho$, $J/\psi$, and $\Upsilon$ mesons were
obtained as outlined in Sect. III.1. The numerical results are

$\mathcal{H}=-ig\begin{cases}84\ \textrm{MeV}^{2}/m_{q},&\rho\\\ 190\
\textrm{MeV}^{2}/{m_{c}},&J/\psi\\\ 225\
\textrm{MeV}^{2}/{m_{b}},&\Upsilon\end{cases}\approx-i\begin{cases}210\
\textrm{MeV},&\rho\\\ 60\ \ \textrm{MeV},&J/\psi\\\ 20\ \
\textrm{MeV},&\Upsilon\end{cases}.$ (28)

Model validation is not simple because phenomenological mixing information is
not available (since hybrid mesons of any type are not firmly established).
Comparison to lattice field theory is difficult because the formalism
automatically produces eigenstates over the $q\bar{q}$ and $q\bar{q}g$ Fock
state sectors. There is, however, one lattice computation (that I am aware of)
that uses the nonrelativistic QCD (NRQCD) formalism. This permits defining
bare Fock states and measuring their overlap, in this case due to the operator
$g\bm{\sigma}\cdot{\bf B}/(2m)$Burch:2003zf . The computation is not easy (the
authors note, “Our charmonium results are plagued with systematic errors which
are not easily quantified”), involving a poorly determined renormalization
constant, and difficulties in scale setting. Nevertheless, the authors
estimate a hybrid component of approximately 2.3% in the $J/\psi$ and 0.4% in
the $\Upsilon$. Approximating these as $|\mathcal{H}/(M_{hyb}-M_{V})|^{2}$ and
using the measured charmonium vector mass splitting of approximately 1150
MeVHadronSpectrum:2012gic , gives

$\displaystyle\mathcal{H}_{\scriptstyle NRQCD}$ $\displaystyle\approx$
$\displaystyle 170\ \textrm{MeV}\ (J/\psi)$ $\displaystyle\approx$
$\displaystyle 70\ \textrm{MeV}\ \ (\Upsilon).$

The results of Eq. 28 are approximately a factor of three below these.
Nevertheless, both computations contain large unquantified uncertainties, and
it is encouraging that they are comparable in size and that the ratio of
results does not follow the naively expected inverse quark mass relationship
for either calculation. Over all, we take these results as evidence in favor
of the utility of the constituent gluon/Coulomb gauge model presented here.
Application of the formalism to the isovector vector mesons will be presented
in the following section.

##  VII Hybrid Phenomenology

We examine the impact of the results presented here on interpreting the light
meson spectrum, with a focus on vector states since these can be made in
$e^{+}e^{-}$ machines. It is unfortunate that experimental knowledge of the
excited rho spectrum is spotty. For example, the $\rho(1450)$ and $\rho(1700)$
are both seen in $4\pi$ or $a_{1}\pi$ decay modes, which is a signal for
hybrid structure. Clarifying the situation with further experimental and
lattice field effort is clearly of interest. Similarly, the $\rho(1900)$
region has conflicting signals and complications from $N\bar{N}$
thresholdWorkman:2022ynf . Finally, the BaBar collaboration has measured
$e^{+}e^{-}\to 2(\pi^{+}\pi^{-})\pi^{0}$, which reveals interesting (although
low statistics) features near 2100 MeVBaBar:2007qju .

Table 4: Model Assignments and Experimental Isovector Vector States (MeV). state | Ref. | mass | width | model (iv) | mass | model (v) | mass
---|---|---|---|---|---|---|---
$\rho(770)$ | PDG | $775.2\pm 0.2$ | $147.4\pm 0.8$ | $1{}^{3}S_{1}$ | 720 | $1{}^{3}S_{1}$ | 810
$\rho(1450)$ | PDG | $1465\pm 25$ | $400\pm 60$ | $2{}^{3}S_{1}$ | 1440 | $2{}^{3}S_{1}$ | 1405
$\rho(1570)$ | PDG | $1570\pm 70$ | $144\pm 90$ | $1{}^{3}D_{1}$ | 1510 | $1{}^{3}D_{1}$ | 1497
$\rho(1700)$ | PDG | $1720\pm 20$ | $250\pm 100$ | $H_{1}(1^{--})$ | 1760 | $3{}^{3}S_{1}$ | 1770
|  |  |  |  |  | $2{}^{3}D_{1}$ | 1835
|  |  |  | $3{}^{3}S_{1}$ | 1850 |  |
$\rho(1900)$ | BaBar:2007ceh  | $1900\pm 30$ | $50\pm 30$ | $2{}^{3}D_{1}$ | 1910 | $4{}^{3}S_{1}$ | 2080
$\rho(2150)$ | BESIII:2020xmw  | $2034\pm 16$ | $234\pm 39$ | $4{}^{3}S_{1}$ | 2170 | $H_{1}(1^{--})$ | 2100
|  |  |  | $3{}^{3}D_{1}$ | 2220 | $3{}^{3}D_{1}$ | 2130

A summary of possible quark model identifications for the rho spectrum is
shown in Table 4. Models (iv) and (v) masses are in rough agreement, but
notice that deviations of tens of MeV low in the spectrum become 90 MeV by the
$4S$ state. This, seemingly minor, difference can lead to substantial changes
in interpretation. For example, there is no vector state near 1700 MeV in
model (iv), raising the possibility that the $\rho(1700)$ may be a hybrid
state. The large symbols in Figure 5 are from a lattice computation at
$m_{\pi}=391$ MeVDudek:2013yja . The $1^{-+}$ state is measured at 2026 MeV,
requiring a shift of 430 MeV to bring it to agreement with a presumed
$\pi_{1}(1600)$. This then implies that the $H_{1}(^{--})$ should have a mass
of 1760 MeV. In view of this is it tempting to make the particle
identifications shown in column 5 of Table 4.

The smaller points in Fig. 5 (obtained from Ref. Dudek:2011bn ) show $H_{1}$
hybrid masses computed at $m_{\pi}=702$, 524, 444, and 396 MeV. These data
permit a rough extrapolation to the physical light quark mass (evidently the
data are not in the chiral regime, so a simple extrapolation in $m_{\pi}$ is
used here). We estimate a $H_{1}(1^{-+})$ mass of 1750 MeV, which implies an
additional shift of 150 MeV. Performing the same procedure for the
$H_{1}(1^{--})$ gives an estimated physical mass of 2100 MeV for the vector
hybrid. A natural candidate for this state is the $\rho(2150)$. The $3S$ mass
in model (v) matches the $\rho(1700)$ well, and we arrive at the alternative
scenario spelled out in the last two columns of Table 4.

There is a clear moral to this story: particle identification relies on
correctly interpreting mass differences of order 100 MeV. Both lattice and
model variations can easily generate deviations of this magnitude, thus it is
important to track model sensitivity in making assignments.

In spite of the lack of moral clarity, it is likely that model (v) is more
robust than model (iv) because it is fit to a much larger array of meson
masses. The model (iv) scenario also has a missing $3S$ state, which is more
unlikely than the missing $2D$ state of model (v). It also seems clear that
the lattice $1^{-+}$ mass has a larger quark mass-dependence than the $1^{--}$
mass. Lastly, as mentioned above, intriguing structure is seen in
$e^{+}e^{-}\to 2(\pi^{+}\pi^{-})\pi^{0}$ near 2100 MeVBaBar:2007qju . In sum,
we feel that the scenario of model (v) should be taken seriously.

Figure 5: Lattice $H_{1}$ MassesDudek:2011bn  (large symbols Ref.
Dudek:2013yja ) vs. Pion Mass. The $\pi_{1}$ is indicated at lower left.

If the $\pi_{1}(1600)$ is confirmed as a hybrid meson and a $\rho(2100)$ is
found that matches expectations for hybrid production and decay, it is natural
to inquire into the accompanying flavor multiplets that must exist. Taken
together, the model and lattice results of Sect. V imply that an isoscalar
$\eta_{1}$ should exist at 1750-1780 MeV, while the “$s\bar{s}$” isoscalar
should have a mass of approximately 1900 MeV. Interestingly, the BESIII
collaboration recently announced the discoveryBESIII:2022riz  of an exotic
isoscalar meson with $J^{PC}=1^{-+}$ quantum numbers, a Breit-Wigner mass of
$1855\pm 9\pm 4$ MeV, and a width of $188\pm 18\pm 5$ MeV. Given the
discussion it is natural to identify this state with the $s\bar{s}$ isoscalar
partner of the $\pi_{1}(1600)$, and suggests that searching for
$\eta_{1}(1760)$ is in order.

As discussed in Sect. VI, the $H_{1}(1^{--})$ hybrid is expected to mix with
nearby canonical vectors. This mixing is explored here to set expectations for
the isovector vector hybrid, whose lattice mass we have suggested is near 2100
MeV. Using typical values for the strong coupling and light quark masses gave
a mixing matrix element of $\mathcal{H}_{1S}\approx 210$ MeV for the
$\rho(770)$ (see Eq. 28). Repeating the calculation for radially excited rho
mesons gives

$\displaystyle\mathcal{H}_{1S}\approx 210\
\textrm{MeV}\qquad\mathcal{H}_{2S}\approx-130\
\textrm{MeV}\qquad\mathcal{H}_{3S}\approx 68\ \textrm{MeV}$
$\displaystyle\mathcal{H}_{4S}\approx-35\
\textrm{MeV}\qquad\mathcal{H}_{5S}\approx 16\
\textrm{MeV}\qquad\mathcal{H}_{6S}\approx-8\ \textrm{MeV}.$

The effect of this mixing on the spectrum is obtained by diagonalizing a
matrix with the diagonals set by the quark model masses of model (v) (see
Table 4). (Strictly speaking bare quantities should appear in the mixing
matrix. But the lattice hybrid mass (should) already account for mixing, while
the quark model should absorb the effects of this mixing (where possible) in
its parameters. However, as will be demonstrated next, mixing is small and
hence this procedure serves as a useful illustration of the expected size of
the effect). The hybrid entry is set to 2100 MeV, while the off-diagonal
entries are set to $\mathcal{H}_{nS}$ in the hybrid:$(nS)$ entry. The
resulting spectrum is displayed versus the strong coupling in Fig. 6. It
appears that mixing effects are small, with the chief outcome being that the
$\rho(770)$ and $H_{1}(1^{--})$ repel as the coupling is increased.
Experimental masses and widths are displayed as boxes and vertical bars in the
figure. The agreement between expectations and experiment is reasonable. Note
especially that the $4S-H_{1}(1^{--})$ splitting is comparable to the
$\rho(2150)-\rho(1900)$ splitting near $g=1.0$, lending support to the “model
(v)” scenario presented above.

With the isovector vector hybrid mass estimated to be 2100 MeV, the lattice
and model results of Sect. V imply that an isoscalar vector hybrid should have
a mass of 2100-2250 MeV, while the “$s\bar{s}$” isoscalar should lie in the
range 2220-2350 MeV.

Figure 6: Vector isovector masses with hybrid mixing as a function of the
strong coupling. Boxes indicate experimental masses and their uncertainties.
Grey bars indicate state widths.

If a vector hybrid is to be discovered in $e^{+}e^{-}$ annihilation its decay
constant should be comparable to other excited $\rho$ states. Presumably this
coupling is set by the $q\bar{q}$ content of the hybrid, which can be obtained
from the mixing matrix just described. The bare $q\bar{q}$ component of the
full hybrid is shown as a function of the strong coupling in Fig. 7, where it
is seen that approximately 20% of the state is $q\bar{q}$ in various
configurations.

Figure 7: $q\bar{q}g$ Fraction of Vector Hybrid.

This observation can be used to compute the hybrid decay constant. Defining a
decay constant in the usual way

$\langle
0|\bar{\psi}\gamma^{\mu}\psi(0)|V(\lambda)\rangle=m_{V}f_{V}\epsilon^{\mu}(\lambda)$
(29)

permits one to obtain the decay (neglecting electron masses)

$\Gamma(V\to e^{+}e^{-})=\frac{4\pi\alpha^{2}Q_{eff}^{2}}{3m_{V}}f_{V}^{2}$
(30)

where $Q_{eff}=Q_{u}/\sqrt{2}+Q_{d}/\sqrt{2}$ is the effective charge of an
isovector meson in units of the electron charge. In the case of the hybrid
meson, allowing the state to be a sum over components gives

$f_{H}=\frac{1}{\sqrt{M_{H}}}\sum_{n\neq H}\sqrt{M_{n}}f_{V}^{(n)}C_{n}$ (31)

where $C_{n}=\langle nS|H_{1}(1^{--})\rangle$ are the state components
obtained previously.

Unfortunately, only the $\rho(770)$ decay constant is known, so we must model
the remaining. Following the methodology of Ref. Lakhina:2006vg , the decay
constant is written as

$f_{V}^{(n)}=\sqrt{\frac{3}{M_{n}}}\int\frac{d^{3}k}{(2\pi)^{3}}\psi^{(n)}(\vec{k})\sqrt{1+\frac{m_{q}}{E_{k}}}\sqrt{1+\frac{m_{\bar{q}}}{E_{\bar{k}}}}\left(1+\frac{k^{2}}{3(E_{k}+m_{q})(E_{\bar{k}}+m_{\bar{q}})}\right).$
(32)

The mass in this expression originates in the relativistic normalization of
the state vector, and is the reason the square roots of meson mass appear in
Eq. 31.

Evaluating Eq. 32 for the $\rho(770)$ gives $f_{\rho}=300$ MeV, to be compared
to the experimental value of 220 MeV (the simple quark model is known to give
decay constants that are too large). Other decay constants obtained in this
way are $f_{\rho(2S)}=160$ MeV, $f_{\rho(3S)}=130$ MeV, $f_{\rho(4S)}=110$
MeV, $f_{\rho(5S)}=100$ MeV, and $f_{\rho(6S)}=95$ MeV. Evaluating Eq. 31 then
yields

$f_{H_{1}(1^{--})}\approx 20\ \textrm{MeV}.$ (33)

where we have accounted for the tendency to over-predict decay constants. Thus
one can expect production of the vector hybrid in $e^{+}e^{-}$ annihilation at
approximately 1% of the strength of canonical mesons.

##  VIII Conclusions

A constituent gluon model of gluodynamics has been explored. This model
permits describing glueballs and hybrid mesons as simple bound states in a
formalism that can be considered a minimal extension of typical constituent
quark models. The model is commensurate with lattice field theory results and
leverages the Hamiltonian of QCD in Coulomb gauge to describe the relevant
dynamics.

This picture was used to model hybrid meson flavor mixing, assuming that
hybrids and glueballs are dominated by the minimal number of quasigluons
required by the state. Hybrid flavor mixing is unique in that it can be
described by low order diagrams because the quark-antiquark pair is in a color
octet state. In this case annihilation can happen via a transverse gluon or an
instantaneous (potential) gluon. Both cases require quark annihilation in the
${}^{(2S+1)}L_{J}={}^{3}S_{1}$ state, which in turn restricts substantial
flavor mixing to the $H_{1}(0^{-+})$, $H_{1}(1^{-+})$, and $H_{1}(2^{-+})$
multiplets. Comparison to lattice results computed with a pion of mass 391 MeV
show broad agreement, with the largest discrepancy in the vector multiplet.
This discrepancy is difficult to reconcile in the present model, and if
confirmed, likely implies that the quasigluon approximation needs to be
abandoned or heavily modified. Of course, any new model must continue to
explain the weak mixing computed in the $H_{2}$ (and presumably other)
multiplets. Alternatively, the model agrees broadly with lattice data,
implying that the assumptions made may be reasonable, and that quasigluons do
indeed serve as a useful description of low-lying gluonic excitations.

A similar computation of configuration mixing of hybrid and canonical mesons
yielded reasonable agreement with an NRQCD computation for $J/\psi$ and
$\Upsilon$ mixing (within large errors for both methods). Thus the quasigluon
approach proves useful in this context as well.

Combining the model computations leads to a picture in which an isovector
vector hybrid is expected with a mass of approximately 2100 MeV. This vector
does not mix extensively with canonical mesons, and has a decay constant of
approximately 20 MeV, perhaps permitting its observation in $e^{+}e^{-}$
scattering. The state’s isospin partners are expected at 2100-2250 MeV and
2220-2350 MeV. Similarly, if the $\pi_{1}(1600)$ is confirmed then isospin
partners are expected at 1750-1780 MeV and around 1900 MeV. This last state is
a natural identification of the recently seen $\eta_{1}$.

The results of this investigation encourage further work on the hybrid
spectrum with the goal of achieving detailed agreement. This will likely
require the addition of higher order spin-dependent and spin-independent
interactions. It will also be of interest to compute light hybrid strong decay
rates as these will be crucial to identifying novel hadrons. Finally, the
topic of radiative transitions of hybrid mesons remains relatively unexplored.

###### Acknowledgements.

The author thanks Matthew Shepherd for rekindling his interest in light
exotics and acknowledges support by the U.S. Department of Energy under
contract DE-SC0019232.

## References

  * (1) C. Meyer and E. Swanson, Prog. Part. Nucl. Phys. 82, 21-58 (2015) doi:10.1016/j.ppnp.2015.03.001 [arXiv:1502.07276 [hep-ph]].
  * (2) J. J. Dudek, Phys. Rev. D 84, 074023 (2011) doi:10.1103/PhysRevD.84.074023 [arXiv:1106.5515 [hep-ph]].
  * (3) D. Horn and J. Mandula, Phys. Rev. D 17, 898 (1978); A. LeYaouanc, L. Oliver, O. Pène, J.-C. Raynal and S. Ono, Z. Phys. C28, 309 (1985); S. Ishida, H. Sawazaki, M. Oda and K. Yamada, Phys. Rev. D47, 179 (1992);Prog. Theor. Phys. 82, 119 (1989); F. Iddir, A. LeYaouanc, L. Oliver, O. Pène, J.-C. Raynal and S. Ono, Phys. Lett. B205, 564 (1988);
  * (4) A. P. Szczepaniak and E. S. Swanson, Phys. Rev. D 65, 025012 (2002) doi:10.1103/PhysRevD.65.025012 [arXiv:hep-ph/0107078 [hep-ph]].
  * (5) E. S. Swanson and A. P. Szczepaniak, Phys. Rev. D 59, 014035 (1999) doi:10.1103/PhysRevD.59.014035 [hep-ph/9804219].
  * (6) P. Guo, A. P. Szczepaniak, G. Galata, A. Vassallo and E. Santopinto, Phys. Rev. D 78, 056003 (2008) doi:10.1103/PhysRevD.78.056003 [arXiv:0807.2721 [hep-ph]].
  * (7) A. P. Szczepaniak and E. S. Swanson, Phys. Lett. B 577, 61-66 (2003) doi:10.1016/j.physletb.2003.10.008 [arXiv:hep-ph/0308268 [hep-ph]].
  * (8) J. J. Dudek et al. [Hadron Spectrum], Phys. Rev. D 88, no.9, 094505 (2013) doi:10.1103/PhysRevD.88.094505 [arXiv:1309.2608 [hep-lat]].
  * (9) A. P. Szczepaniak and P. Krupinski, Phys. Rev. D 73, 116002 (2006) doi:10.1103/PhysRevD.73.116002 [arXiv:hep-ph/0604098 [hep-ph]].
  * (10) P. Guo, A. P. Szczepaniak, G. Galata, A. Vassallo and E. Santopinto, Phys. Rev. D 77, 056005 (2008) doi:10.1103/PhysRevD.77.056005 [arXiv:0707.3156 [hep-ph]].
  * (11) C. Farina, H. Garcia Tecocoatzi, A. Giachino, E. Santopinto and E. S. Swanson, Phys. Rev. D 102, no.1, 014023 (2020) doi:10.1103/PhysRevD.102.014023 [arXiv:2005.10850 [hep-ph]].
  * (12) N. H. Christ and T. D. Lee, Phys. Rev. D 22, 939 (1980).
  * (13) C. Feuchter and H. Reinhardt, Phys. Rev. D 70, 105021 (2004) doi:10.1103/PhysRevD.70.105021 [arXiv:hep-th/0408236 [hep-th]].
  * (14) Y. Koma and M. Koma, “Spin-dependent potentials from lattice QCD,” Nucl. Phys. B 769, 79 (2007) [hep-lat/0609078].
  * (15) K. Juge, J. Kuti and C. Morningstar, Phys. Rev. Lett. 82, 4400-4403 (1999) doi:10.1103/PhysRevLett.82.4400 [arXiv:hep-ph/9902336 [hep-ph]].
  * (16) See Chapter 21 in J.F. Boudreau and E.S. Swanson, Applied Computational Physics, (Oxford U. P., Oxford, 2018).
  * (17) R. L. Workman et al. [Particle Data Group], PTEP 2022, 083C01 (2022) doi:10.1093/ptep/ptac097
  * (18) A. Rodas et al. [JPAC], Phys. Rev. Lett. 122, no.4, 042002 (2019) doi:10.1103/PhysRevLett.122.042002 [arXiv:1810.04171 [hep-ph]].
  * (19) B. Kopf, M. Albrecht, H. Koch, M. Küßner, J. Pychy, X. Qin and U. Wiedner, Eur. Phys. J. C 81, no.12, 1056 (2021) doi:10.1140/epjc/s10052-021-09821-2 [arXiv:2008.11566 [hep-ph]].
  * (20) L. Liu et al. [Hadron Spectrum], JHEP 07, 126 (2012) doi:10.1007/JHEP07(2012)126 [arXiv:1204.5425 [hep-ph]].
  * (21) G. K. C. Cheung et al. [Hadron Spectrum], JHEP 12, 089 (2016) doi:10.1007/JHEP12(2016)089 [arXiv:1610.01073 [hep-lat]].
  * (22) H.J. Lipkin, Phys Lett B124, 509 (1983) doi:10.1016/0370-2693(83)91563-0
  * (23) P. Geiger and N. Isgur, Phys. Rev. D 44, 799-808 (1991) doi:10.1103/PhysRevD.44.799
  * (24) P. Geiger and N. Isgur, Phys. Rev. D 47, 5050-5059 (1993) doi:10.1103/PhysRevD.47.5050
  * (25) O. Lakhina and E. S. Swanson, Phys. Rev. D 74, 014012 (2006) doi:10.1103/PhysRevD.74.014012 [arXiv:hep-ph/0603164 [hep-ph]].
  * (26) S. J. Brodsky, C. R. Ji, A. Pang and D. G. Robertson, Phys. Rev. D 57, 245-252 (1998) doi:10.1103/PhysRevD.57.245 [arXiv:hep-ph/9705221 [hep-ph]].
  * (27) Y. Chen, A. Alexandru, S. J. Dong, T. Draper, I. Horvath, F. X. Lee, K. F. Liu, N. Mathur, C. Morningstar and M. Peardon, et al. Phys. Rev. D 73, 014516 (2006) doi:10.1103/PhysRevD.73.014516 [arXiv:hep-lat/0510074 [hep-lat]].
  * (28) T. Barnes, Nucl. Phys. B 158, 171-188 (1979) doi:10.1016/0550-3213(79)90194-9
  * (29) S. B. Gerasimov, [arXiv:hep-ph/9812509 [hep-ph]].
  * (30) T. Burch et al. [MILC], Phys. Rev. D 68, 094504 (2003) doi:10.1103/PhysRevD.68.094504 [arXiv:hep-lat/0305008 [hep-lat]].
  * (31) B. Aubert et al. [BaBar], Phys. Rev. D 76, 092005 (2007) [erratum: Phys. Rev. D 77, 119902 (2008)] doi:10.1103/PhysRevD.76.092005 [arXiv:0708.2461 [hep-ex]].
  * (32) B. Aubert et al. [BaBar], Phys. Rev. D 77, 092002 (2008) doi:10.1103/PhysRevD.77.092002 [arXiv:0710.4451 [hep-ex]].
  * (33) M. Ablikim et al. [BESIII], Phys. Lett. B 813, 136059 (2021) doi:10.1016/j.physletb.2020.136059 [arXiv:2009.08099 [hep-ex]].
  * (34) M. Ablikim et al. [BESIII], Phys. Rev. Lett. 129, no.19, 192002 (2022) doi:10.1103/PhysRevLett.129.192002 [arXiv:2202.00621 [hep-ex]].


# Integrable maps in 4D and modified Volterra lattices

A.N.W. Hone and J.A.G. Roberts
School of Mathematics & Statistics
University of New South Wales
Sydney, NSW 2052, Australia.
<EMAIL_ADDRESS><EMAIL_ADDRESS>

P. Vanhaecke
Laboratoire de Mathématiques et Applications
UMR 7348 CNRS-Université de Poitiers
86360 Chasseneuil-du-Poitou, France.
<EMAIL_ADDRESS>

F. Zullo
DICATAM, Università degli Studi di Brescia
via Branze, 38 - 25123 Brescia, Italia &
INFN, sezione di Milano-Bicocca
Piazza della Scienza, 3 20126, Milano, Italy
<EMAIL_ADDRESS>Work begun while on leave from School of Mathematics,
Statistics & Actuarial Science, University of Kent, Canterbury CT2 7NF, UK.

###### Abstract

In recent work, we presented the construction of a family of difference
equations associated with the Stieltjes continued fraction expansion of a
certain function on a hyperelliptic curve of genus $g$. As well as proving
that each such discrete system is an integrable map in the Liouville sense, we
also showed it to be an algebraic completely integrable system. In the
discrete setting, the latter means that the generic level set of the
invariants is an affine part of an abelian variety, in this case the Jacobian
of the hyperelliptic curve, and each iteration of the map corresponds to a
translation by a fixed vector on the Jacobian. In addition, we demonstrated
that, by combining the discrete integrable dynamics with the flow of one of
the commuting Hamiltonian vector fields, these maps provide genus $g$ algebro-
geometric solutions of the infinite Volterra lattice, which justified naming
them Volterra maps, denoted ${\cal V}_{g}$.

The original motivation behind our work was the fact that, in the particular
case $g=2$, we could recover an example of an integrable symplectic map in
four dimensions found by Gubbiotti, Joshi, Tran and Viallet, who classified
birational maps in 4D admitting two invariants (first integrals) with a
particular degree structure, by considering recurrences of fourth order with a
certain symmetry. Hence, in this particular case, the map ${\cal V}_{2}$
yields genus two solutions of the Volterra lattice. The purpose of this note
is to point out how two of the other 4D integrable maps obtained in the
classification of Gubbiotti et al. correspond to genus two solutions of two
different forms of the modified Volterra lattice, being related via a Miura-
type transformation to the $g=2$ Volterra map ${\cal V}_{2}$.

We dedicate this work to a dear friend and colleague, Decio Levi.

##  1 Introduction

This short article consists of some recollections of our colleague Decio Levi
(in section 2 below), followed by a brief update on our recent results about
integrable maps in four (and higher) dimensions, which provide algebro-
geometric solutions of differential-difference equations of Volterra type
[10]. Decio was one of the pioneers in the theory of integrability for
differential-difference equations, especially in the construction of
integrable lattices from Bäcklund transformations for continuous systems [13,
15], and the programme of applying the symmetry approach to the classification
of such lattices, which he initiated with Yamilov [14]. Thus we like to think
that Decio would have appreciated the results being presented here.

After presenting a few memories of Decio, in section 3 we begin by giving a
brief overview of the 4D integrable maps which were classified by Gubbiotti et
al. [7]. We then proceed to review our construction of integrable maps
obtained from the Stieltjes fraction expansion of certain functions on
hyperelliptic curves [10], and explain how it reproduces one of the examples
from [7], denoted (P.iv), in the particular case of genus two curves. Sections
4 and 5 are devoted to the maps (P.v) and (P.vi), respectively: we show how
each of these maps is related to a different form of the modified Volterra
lattice, and present explicit formulae which relate their solutions to the
solutions of (P.iv) via a transformation of Miura type. We end with some very
short conclusions in section 6.

##  2 Memories of Decio Levi

Andrew Hone writes:  I first met Decio in Warsaw in September 1995, when I was
a PhD student participating in the 1st Non-Orthodox School on Nonlinearity and
Geometry [23]. Decio was one of the lecturers, along with Orlando Ragnisco,
and it was thanks to extended conversations with Orlando that I resolved to
apply for postdoctoral funding to work with him when I finished my PhD. After
receiving a grant from the Leverhulme Trust two years later, I finally got to
be a researcher at Roma Tre, where Orlando and Decio were both professors in
the Dipartimento di Fisica.

For approximately the first six months of my time in Rome, there was no
available office space for postdocs, which meant that I had to share an office
with Orlando. Far from being a negative aspect of my experience, this
situation had many positive benefits for me, and not just scientific ones. By
working in close proximity with Orlando, it meant that I was privy to the
regular visits from the neighbour in the office next door, namely Decio, his
long-time friend and collaborator. Apart from the pleasure of getting to know
Decio, and learning many wonderful ideas about integrable systems from him,
there was the fact that, by default, he would chat to Orlando in Italian,
which helped me to rapidly improve my grasp of the language in those first few
months. The strong bond of friendship between Orlando and Decio created a very
happy atmosphere, and I have extremely fond memories of those times.

In subsequent years, I would see Decio fairly often at various international
conferences, or during return visits to Rome. He had an amiable manner and a
warm, cheerful smile. It was always enjoyable to talk to him, whether about
technical problems, sharing family news, or just musing about life in general.
Talking with Decio would leave me feeling reassured, that all was right with
the world, and I liked his gentle way of concluding a long conversation with
“Vabbè in somma”.

It is an honour to be able to remember Decio here, both for his contributions
as a scientist, and as a wonderful human being.

Federico Zullo writes:  The first time I met Decio was in 2003: I was a
student at the Dipartimento di Fisica of Roma Tre University and needed an
advisor for my last examination for my laurea triennale (bachelor’s degree). I
asked Orlando Ragnisco who, at that time, was very busy. He accompanied me to
the office next door, where Decio was, and I asked him for a theme for my
short dissertation. He very heartily introduced me to the subject of solitons,
that I never heard about before, giving me books and kind advice. Later,
during my laurea magistrale (master’s degree), and during my PhD studies, I
followed different classes taught by Decio, some with very few students. The
familiar atmosphere and natural mildness of Decio’s classes fostered my
learning, and I’m greatly indebted to him for having taught me many topics
used in mathematical physics, like group theory, symmetries of differential
equations, physics of nonlinear systems, qualitative and quantitative analysis
of solutions of differential equations and others. For my own teaching, I
still use some of the material that I collected from his courses.

For a period just before 2014, I was hosted by Decio in his office as a
researcher. I remember the talks on disparate subjects, like religion,
literature, politics, society and, obviously, our research. The talks would
then continue during the lunch break, usually in Via Marconi, with Orlando and
the other members of the very stimulating group of young researchers that was
gathered at Roma Tre in that period, including Fabio Musso, Matteo Petrera,
Christian Scimiterna, Danilo Riglioni, Riccardo Droghei, and later Giorgio
Gubbiotti and Danilo Latini, all led by Decio and Orlando. I’ll always keep
these beautiful memories with me.

##  3 The map (P.iv) and the geometry of its solutions

Discrete integrable systems can be constructed by applying an appropriate
discretization procedure to continuous ones, and historically this is how many
examples of discrete integrability were first discovered [13, 20]. However,
from both a theoretical point of view and a practical one, it is important to
have a notion of integrability for discrete systems that does not require
making reference to some underlying continuous system, whether this be for
lattice equations [14], or for integrable maps [2, 16, 22]. While integrable
maps in two and three dimensions lead to families of invariant curves (as the
level sets of first integrals), the case of four dimensions can lead to new
features, namely invariant tori of dimension two.

In [7], Gubbiotti et al. presented a classification of four-dimensional
birational maps of recurrence type, that is

$\varphi:\qquad(w_{0},w_{1},w_{2},w_{3})\mapsto\Big{(}w_{1},w_{2},w_{3},F(w_{0},w_{1},w_{2},w_{3})\Big{)},$
(3.1)

for a suitable rational function $F$ of the affine coordinates
$(w_{0},w_{1},w_{2},w_{3})\in{\mathbb{C}}^{4}$, where the map $\varphi$ is
required to be invariant under the involution
$\iota:\,(w_{0},w_{1},w_{2},w_{3})\mapsto(w_{3},w_{2},w_{1},w_{0})$, and to
possess two independent polynomial invariants, $H_{1}$, $H_{2}$ say, with
specific degree patterns
$(\deg_{w_{0}}H_{j},\deg_{w_{1}}H_{j},\deg_{w_{2}}H_{j},\deg_{w_{3}}H_{j})=(1,3,3,1)$
and $(2,4,4,2)$ for $j=1,2$, respectively. The result of this classification
was six maps with parameters, labelled (P.i-vi), together with six associated
maps, denoted (Q.i-vi) respectively. Each of the “Q” maps arises from a
corresponding “P” map, as a discrete integrating factor for linear
combinations of the first integrals, so they are dual to one another in the
sense of [17].

As described previously, first in [11] and then [6], the original motivation
for classifying such maps was to understand autonomous versions of the fourth-
order members of hierarchies of discrete Painlevé I/II equations from [5];
but, aside from the latter connection, the “P” in this nomenclature has
nothing to do with the usual labelling of continuous Painlevé equations. From
our point of view, the most interesting cases are the maps labelled (P.iv),
(P.v) and (P.vi), since (from Table 1 in [7]) these are the only ones arising
from a discrete variational principle (Lagrangian), leading to a non-
degenerate Poisson bracket in four dimensions, such that the two first
integrals $H_{1}$, $H_{2}$ are in involution; this means that in the real case
the Liouville tori are two-dimensional. Subsequently, Gubbiotti obtained these
4D integrable maps via an alternative method, by classifying fourth-order
difference equations with a discrete Lagrangian structure [8].

Here we begin with the case of (P.iv), which is the birational map given in
affine coordinates by the recurrence

$\begin{array}[]{l}w_{n+4}w_{n+3}w_{n+2}+w_{n+2}w_{n+1}w_{n}+2w_{n+2}^{2}(w_{n+3}+w_{n+1})\\\
+w_{n+2}(w_{n+3}^{2}+w_{n+3}w_{n+1}+w_{n+1}^{2})+w_{n+2}^{3}+\nu
w_{n+2}(w_{n+3}+w_{n+2}+w_{n+1})+bw_{n+2}+a=0.\end{array}$ (3.2)

This map has three essential parameters $a,b,\nu$ (in the formulae from [7] we
have set the parameter $d=1$, which can be achieved by a simple rescaling),
and it is of the form (3.1), with

$F=-\frac{w_{0}w_{1}w_{2}+w_{1}w_{2}w_{3}+w_{1}^{2}w_{2}+w_{2}w_{3}^{2}+2w_{1}w_{2}^{2}+2w_{2}^{2}w_{3}+w_{2}^{3}+\nu(w_{1}w_{2}+w_{2}w_{3}+w_{2}^{2})+bw_{2}+a}{w_{2}w_{3}};$

this $F$ is the rational function of $w_{0},w_{1},w_{2},w_{3}$ obtained by
solving for $w_{4}$ in (3.2) with $n=0$.

The first integral denoted $I_{\mathrm{low}}^{\mathrm{P}.\mathrm{iv}}$ in [7]
is given in affine coordinates by

$H_{1}=w_{1}w_{2}\Big{(}w_{2}w_{3}+w_{0}w_{1}-w_{0}w_{3}+(w_{1}+w_{2})^{2}+\nu(w_{1}+w_{2})+b\Big{)}+a(w_{1}+w_{2}).$
(3.3)

The latter has the degree pattern $(1,3,3,1)$. In particular, it is linear in
$w_{3}$, which implies that, on each three-dimensional level set
$H_{1}=h_{1}=\,\,$const, the map (3.2) reduces to a birational map in three
dimensions, given by the recurrence

$\begin{array}[]{rcl}w_{n+3}w_{n+2}w_{n+1}(w_{n+2}-w_{n})+w_{n+2}w_{n+1}^{2}w_{n}+w_{n+2}w_{n+1}(w_{n+1}+w_{n+2})^{2}&&\\\
+\nu\,w_{n+2}w_{n+1}(w_{n+1}+w_{n+2})+b\,w_{n+2}w_{n+1}+a\,(w_{n+1}+w_{n+2})&=&h_{1}.\end{array}$

A second independent invariant for (3.2), with degree pattern $(2,4,4,2)$, is
given by

$\begin{array}[]{rcl}H_{2}&=&w_{1}w_{2}\left(\begin{array}[]{c}w_{0}^{2}w_{1}+w_{3}^{2}w_{2}+w_{0}w_{3}(w_{1}+w_{2})+w_{0}(w_{2}^{2}+2w_{1}^{2})+w_{3}(w_{1}^{2}+2w_{2}^{2})\\\
+\,3(w_{0}+w_{3})w_{1}w_{2}+(w_{1}+w_{2})^{3}\\\
+\nu\,\,\Big{(}w_{0}w_{3}+(w_{0}+w_{3})(w_{1}+w_{2})+(w_{1}+w_{2})^{2}\Big{)}+b\,(w_{0}+w_{1}+w_{2}+w_{3})\end{array}\right)\\\
&&+a\,\Big{(}w_{0}w_{1}+w_{3}w_{2}+(w_{1}+w_{2})^{2}\Big{)}.\end{array}$ (3.4)

This differs slightly from the second invariant presented in [7], which is
$I_{\mathrm{high}}^{\mathrm{P}.\mathrm{iv}}=H_{2}-\nu H_{1}$.

The nondegenerate Poisson bracket between the coordinates, which was obtained
in [7] by making use of a discrete Lagrangian for (3.2), is given by

$\\{\,w_{n},w_{n+1}\,\\}=0,\,\\{\,w_{n},w_{n+2}\,\\}=\frac{1}{w_{n+1}},\,\\{\,w_{n},w_{n+3}\,\\}=-\frac{w_{n}+2w_{n+1}+2w_{n+2}+w_{n+3}+\nu}{w_{n+1}w_{n+2}},$
(3.5)

for all $n$. So (3.1) is a Poisson map, in the sense that
$\\{\,\varphi^{*}G,\varphi^{*}H\,\\}=\varphi^{*}\\{\,G,H\,\\}$ for all
functions $G,H$ on ${\mathbb{C}}^{4}$. The two independent invariants given in
[7] are in involution with respect to this bracket, which is equivalent to the
involutivity of functions (3.3) and (3.4), that is to say

$\\{\,H_{1},H_{2}\,\\}=0.$

Hence the four-dimensional map defined by (3.2) is integrable in the Liouville
sense.

Computing the Hamiltonian vector field for the first flow, generated by
$H_{1}$, we find that this takes the form

$\frac{\mathrm{d}w_{n}}{\mathrm{d}t}=w_{n}(w_{n+1}-w_{n-1})$ (3.6)

for $n=1,2$. However, since (3.2) is a Poisson map that commutes with this
flow, it follows that the relation (3.6) extends to all $n\in{\mathbb{Z}}$.
Thus the combined solutions of the map and the flow, which are compatible with
one another, generate a sequence of functions
$\big{(}w_{n}(t)\big{)}_{n\in{\mathbb{Z}}}$ satisfying (3.6), which is the
Volterra lattice equation, first considered by Kac and van Moerbeke [12].
Hence, in a certain sense that can be made precise, these will turn out to be
genus 2 solutions of this lattice hierarchy.

The complex geometry of the solutions of the map defined by (3.2) is related
to a family of hyperelliptic curves of genus 2, given by the Weierstrass
quintic

$\Gamma:\quad y^{2}=(1+\nu x+bx^{2})^{2}+4a(1+\nu
x+bx^{2})x^{3}+4h_{1}x^{4}+4(h_{2}+\nu h_{1})x^{5}.$ (3.7)

On any genus 2 curve $\Gamma$ of the above form, we take the meromorphic
function $F$ given by

$F=\frac{y+{\cal P}(x)}{{\cal Q}(x)}=\frac{{\cal R}(x)}{y-{\cal P}(x)},$ (3.8)

where ${\cal P},{\cal Q},{\cal R}$ are polynomials in the spectral parameter
$x$, given by

${\cal P}(x)=1+p_{1}x+p_{2}x^{2},\quad{\cal
Q}(x)=2+q_{1}x+q_{2}x^{2},\quad{\cal R}(x)=r_{1}x+r_{2}x^{2}+r_{3}x^{3},$
(3.9)

which are required to satisfy

${\cal P}(x)^{2}+{\cal Q}(x){\cal R}(x)=f(x),$ (3.10)

with $f(x)=(1+\nu x+bx^{2})^{2}+4a(1+\nu
x+bx^{2})x^{3}+4h_{1}x^{4}+4(h_{2}+\nu h_{1})x^{5}$ being the quintic on the
right-hand side of (3.7). Then the key to the construction in [10] is to
expand the function $F$ as a continued fraction of Stieltjes type
(S-fraction), that is

$F=1-\cfrac{w_{1}x}{1-\cfrac{w_{2}x}{1-\cfrac{w_{3}x}{1-\cdots}}}\;,$ (3.11)

and by iterating from one line of the fraction to the next we find that we
obtain a recurrence for the coefficients $w_{j}$. More precisely, the non-
trivial coefficients of the polynomials (3.9) are given in terms of $w_{j}$
and the parameters by

$p_{1}=2w_{0}+\nu,\,p_{2}=2w_{0}(w_{1}+w_{0}+w_{-1})+b,\,\tfrac{1}{2}q_{1}=w_{0}+w_{1}+\nu,\,\tfrac{1}{2}q_{2}=w_{0}w_{-1}+w_{1}w_{2}+(w_{1}+w_{0})^{2}+\nu(w_{0}+w_{1})+b,$

and $r_{1}=-2w_{0}$; there are similar (but slightly more unwieldy)
expressions for $r_{2}$ and $r_{3}$, which are omitted here, but are easily
obtained from the relation (3.10). With these identifications, the iteration
of the S-fraction (3.11) for $F$ becomes precisely the map (P.iv) in terms of
the affine coordinates $w_{j}$, as given by (3.2).

In [10] it was also shown that each iteration of the continued fraction is
equivalent to the discrete Lax equation

${\bf L}(x){\bf M}(x)={\bf M}(x)\widetilde{\bf L}(x)\;,$

where

${\bf L}(x):=\left(\begin{array}[]{lr}{\cal P}(x)&{\cal R}(x)\\\ {\cal
Q}(x)&-{\cal P}(x)\end{array}\right)\;,\qquad{\bf
M}(x):=\left(\begin{array}[]{cc}1&-w_{1}x\\\ 1&0\end{array}\right).$

Furthermore, we found that each generic common level set of the two invariants
$H_{1},H_{2}$ is isomorphic to an affine part of the Jacobian of the
associated spectral curve $\Gamma$ (or rather, of its completion), and each
iteration of the map corresponds to a translation on the Jacobian by the
divisor class $[(0,-1)-\infty]$. Thus, in addition to being integrable in the
Liouville sense, the map (3.2) is an algebraic completely integrable system,
being a discrete analogue of an a.c.i. system (see [1, 21]).

The map (3.2) can also be rewritten in terms of tau functions $\tau_{n}$,
related to $w_{n}$ via

$w_{n}=\frac{\tau_{n}\tau_{n+3}}{\tau_{n+1}\tau_{n+2}}.$

These tau functions satisfy a Somos-9 recurrence, that is

${\alpha}_{1}\,\tau_{n+9}\tau_{n}+{\alpha}_{2}\,\tau_{n+8}\tau_{n+1}+{\alpha}_{3}\,\tau_{n+7}\tau_{n+2}+{\alpha}_{4}\,\tau_{n+6}\tau_{n+3}+{\alpha}_{5}\,\tau_{n+5}\tau_{n+4}=0,$
(3.12)

with coefficients ${\alpha}_{j}$ that depend on $a,b,\nu$ and the values of
$H_{1},H_{2}$ along each orbit of (3.2); for details see Proposition 2.1 in
[10]. Using the S-fraction (3.11), we were also able to write explicit Hankel
determinant formulae for these tau functions $\tau_{n}$, analogous to results
for Somos sequences in genus 1 [3], and other Hankel determinant formulae for
solutions of the Volterra lattice [4]. Furthermore, we found a Miura map
relating the solutions of (P.iv) to one of the maps derived from J-fractions
in [9], using the classical method of contraction of continued fractions due
to Stieltjes [19] (see also [18]), which in this case turned out to provide
solutions of the infinite Toda lattice.

In what follows, we will present analogous properties for the maps (P.v) and
(P.vi), and point out how they are closely connected to (P.iv).

##  4 The map (P.v)

The map (P.v) is given by the recurrence

$\begin{array}[]{l}w_{n+4}w_{n+3}^{2}w_{n+2}^{2}+w_{n+2}^{2}w_{n+1}^{2}w_{n}+w_{n+2}^{3}(w_{n+1}+w_{n+3})^{2}+\tilde{\nu}w_{n+2}^{2}(w_{n+1}+w_{n+3})+\tilde{c}w_{n+2}+\tilde{a}=0,\end{array}$
(4.1)

with three essential parameters $\tilde{a},\tilde{c},\tilde{\nu}$ (compared
with [7] we have put tildes here to distinguish them from the parameters in
(3.2), and rescaled so that the parameter $d\to 1$).

The lowest degree first integral of the map defined by (4.1), with degree
pattern $(1,3,3,1)$, is given by

$H_{1}=w_{3}w_{2}^{3}w_{1}^{2}+w_{2}^{2}w_{1}^{3}w_{0}-w_{3}w_{2}^{2}w_{1}^{2}w_{0}+w_{2}^{3}w_{1}^{3}+\tilde{\nu}w_{2}^{2}w_{1}^{2}+\tilde{c}w_{2}w_{1}+\tilde{a}(w_{2}+w_{1}),$
(4.2)

and this is the same as $I_{\mathrm{low}}^{\mathrm{P}.\mathrm{v}}$ in [7].
Another first integral, with degree pattern $(2,4,4,2)$, is

$\begin{array}[]{rcl}H_{2}&=&w_{2}^{2}w_{1}^{2}\Big{(}(w_{3}w_{2}+w_{1}w_{0}+w_{2}w_{1})^{2}+\tilde{\nu}(w_{3}+w_{1})(w_{2}+w_{0})\Big{)}\\\
&&+\tilde{c}w_{2}w_{1}(w_{3}w_{2}+w_{1}w_{0}+w_{2}w_{1})+\tilde{a}(w_{3}w_{2}^{2}+w_{1}^{2}w_{0}+w_{2}^{2}w_{1}+w_{2}w_{1}^{2}).\end{array}$
(4.3)

The second invariant presented in [7] is
$I_{\mathrm{high}}^{\mathrm{P}.\mathrm{v}}=H_{2}-\tilde{\nu}H_{1}$.

The nondegenerate Poisson bracket between the coordinates is given by

$\\{\,w_{n},w_{n+1}\,\\}=0,\,\\{\,w_{n},w_{n+2}\,\\}=\frac{1}{w_{n+1}^{2}},\,\\{\,w_{n},w_{n+3}\,\\}=-\frac{2(w_{n}w_{n+1}+w_{n+1}w_{n+2}+w_{n+2}w_{n+3})+\tilde{\nu}}{w_{n+1}^{2}w_{n+2}^{2}}.$

The independent first integrals (4.2) and (4.3) are in involution with respect
to this bracket, which shows that the map (4.1) is Liouville integrable.

Computing the Hamiltonian vector field for the first flow, generated by
$H_{1}$, we find that this takes the form

$\frac{\mathrm{d}w_{n}}{\mathrm{d}t}=w_{n}^{2}(w_{n+1}-w_{n-1})$ (4.4)

for $n=1,2$. However, since the map (4.1) is Poisson and commutes with the
flow $\\{\cdot,H_{1}\\}$, the equation (4.4) holds for all $n\in{\mathbb{Z}}$.
Thus the compatible solutions of the map and the flow together provide a
sequence of functions $\big{(}w_{n}(t)\big{)}_{n\in{\mathbb{Z}}}$ which
satisfy (4.4), which is a degenerate case of the modified Volterra lattice
equation [25].

If we make the tau function substitution

$w_{n}=\frac{\tau_{n}\tau_{n+2}}{\tau_{n+1}^{2}}$ (4.5)

for (P.v), then we find that the sequence $(\tau_{n})$ satisfies a Somos-8
relation. More precisely, by direct computer algebra calculations we can show
the following:

######  Proposition 4.1.

Whenever $w_{n}$ is a solution of (4.1), the sequence $(\tau_{n})$ satisfies
the following Somos-8 recurrence, with coefficients that are functions of the
Hamiltonians $H_{1},H_{2}$ as in (4.2) and (4.3) above (constant along each
orbit):

${\alpha}_{1}\,\tau_{n+8}\tau_{n}+{\alpha}_{2}\,\tau_{n+7}\tau_{n+1}+{\alpha}_{3}\,\tau_{n+6}\tau_{n+2}+{\alpha}_{4}\,\tau_{n+5}\tau_{3}+{\alpha}_{5}\,\tau_{n+4}^{2}=0,$
(4.6)

where the coefficients are given by

$\begin{array}[]{c}{\alpha}_{1}=H_{1},\qquad{\alpha}_{2}=\tilde{a}H_{2},\qquad{\alpha}_{3}=\tilde{a}^{2}H_{2}-H_{1}^{3},\\\
{\alpha}_{4}=\tilde{a}\Big{(}H_{2}^{2}+\tilde{\nu}H_{1}H_{2}+\tilde{c}H_{1}^{2}+\tilde{a}^{2}H_{1}\Big{)},\qquad{\alpha}_{5}=-H_{1}\Big{(}H_{2}^{2}+\tilde{\nu}H_{1}H_{2}+\tilde{c}H_{1}^{2}+\tilde{a}^{2}H_{1}\Big{)}.\end{array}$

Let us denote a solution of the Volterra lattice (3.6) by $\hat{w}_{n}$. Then
the Miura map from the modified Volterra lattice (4.4) takes the form

$\hat{w}_{n}=w_{n+1}w_{n}.$ (4.7)

This Miura map remains valid at the level of the maps (3.2) and (4.1), in the
following sense.

######  Theorem 4.2.

Let $w_{n}$ be a solution of (4.1) with parameters
$\tilde{a},\tilde{c},\tilde{\nu}$, lying on the level set
$H_{1}=\tilde{h}_{1}$, $H_{2}=\tilde{h}_{2}$, of the first integrals (4.2) and
(4.3). Then $\hat{w}_{n}$ given by the Miura map (4.7) is a solution of (3.2)
with parameters

$\nu=\tilde{\nu},\qquad b=\tilde{c},\qquad a=\tilde{h}_{1}.$

Furthermore, on this solution $\hat{w}_{n}$, the values $h_{1},h_{2}$ of the
first integrals (3.3) and (3.4) for the map (3.2) are given by

$h_{1}=\tilde{h}_{2},\qquad
h_{2}=-\tilde{a}^{2}-\tilde{\nu}\tilde{h}_{2}-\tilde{c}\tilde{h}_{1}.$

  * Proof:

The first part of this result is verified by substituting the Miura formula
(4.7) directly into (3.2), using (4.1) to eliminate $w_{n+5}$ followed by
$w_{n+4}$, and then using the formula for $H_{1}$ in (4.2) to eliminate
$w_{n+3}$ on the level set $H_{1}=\tilde{h}_{1}$. Analogous calculations,
rewriting (3.3) and (3.4) in terms of $w_{n}$ satisfying (4.1) and comparing
with $\tilde{h}_{2}$, the value of the first integral (4.3) for the latter
map, yield the above expressions for $h_{1},h_{2}$. ∎

It is worth commenting on the meaning of the Miura formula (4.7), restricted
to this finite-dimensional setting. Given initial data
$w_{0},w_{1},w_{2},w_{3}$ for the map (4.1), we can fix a level set
$H_{1}=\tilde{h}_{1}$ to write

$\hat{w}_{0}=w_{0}w_{1},\,\hat{w}_{1}=w_{1}w_{2},\,\hat{w}_{2}=w_{2}w_{3},\,\hat{w}_{3}=w_{3}\,G(w_{1},w_{2},w_{3},\tilde{h}_{1}),$

for some rational function $G$, obtained by using the formula (4.2) for
$H_{1}$ to eliminate $w_{4}$. Similarly, we can use $H_{1}$ to eliminate
$w_{0}$ above in terms of $w_{1},w_{2},w_{3}$ and $\tilde{h}_{1}$, and after
taking resultants we can do further elimination to solve for each of
$w_{0},w_{1},w_{2},w_{3}$ as algebraic functions of
$\hat{w}_{0},\hat{w}_{1},\hat{w}_{2},\hat{w}_{3}$ and $\tilde{h}_{1}$. So this
leads to an explicit inverse of (4.7), at least in the form of an algebraic
correspondence.

##  5 The map (P.vi)

The map (P.vi) is given by

$\begin{array}[]{l}w_{n+4}(w_{n+3}^{2}-{\delta}^{2})(w_{n+2}^{2}-{\delta}^{2})+w_{n}(w_{n+1}^{2}-{\delta}^{2})(w_{n+2}^{2}-{\delta}^{2})\\\
+w_{n+2}\Big{(}(w_{n+2}^{2}-{\delta}^{2})(w_{n+3}+w_{n+1})^{2}+\bar{c}-{\delta}^{4}\Big{)}+\bar{\nu}(w_{n+2}^{2}-{\delta}^{2})(w_{n+3}+w_{n+1})+\bar{a}=0.\end{array}$
(5.1)

This depends on only three essential parameters $\bar{a},\bar{c},\bar{\nu}$;
compared with [7] we have replaced $a\to\bar{a}$, $c\to\bar{c}$,
$d\to-\bar{\nu}$ and $\delta\to\delta^{2}$. Note the map P(v) in the previous
section arises from P(vi) in the limit $\delta\to 0$, while for $\delta\neq 0$
the map can always be rescaled so that $\delta\to 1$, but it will be
convenient to retain this parameter which has the same weight as $w_{n}$ in
(5.1).

The lowest degree first integral of the map defined by (5.1), with degree
pattern $(1,3,3,1)$, is given by

$\begin{array}[]{rcl}H_{1}&=&\big{(}w_{1}^{2}w_{2}^{2}-{\delta}^{2}(w_{1}^{2}+w_{2}^{2})\big{)}\Big{(}w_{3}w_{2}+w_{0}w_{1}+w_{1}w_{2}-w_{3}w_{0}+\bar{\nu}\Big{)}\\\
&&{\delta}^{4}(w_{3}w_{2}+w_{0}w_{1}-w_{0}w_{3})+\bar{c}w_{2}w_{1}+\bar{a}(w_{2}+w_{1}).\end{array}$
(5.2)

A nondegenerate Poisson bracket between the coordinates is given by

$\\{\,w_{n},w_{n+1}\,\\}=0,\,\\{\,w_{n},w_{n+2}\,\\}=\frac{1}{w_{n+1}^{2}-{\delta}^{2}},\,\\{\,w_{n},w_{n+3}\,\\}=-\frac{2(w_{n}w_{n+1}+w_{n+1}w_{n+2}+w_{n+2}w_{n+3})+\bar{\nu}}{(w_{n+1}^{2}-{\delta}^{2})(w_{n+2}^{2}-{\delta}^{2})},$

and was derived in [7] using a discrete Lagrangian structure for (5.1). A
second independent first integral $H_{2}$ was given in [7], which is in
involution with $H_{1}$ with respect to this bracket. Here we take the second
independent quantity as

$\begin{array}[]{rcl}H_{2}&=&(w_{1}^{2}-{\delta}^{2})(w_{2}^{2}-{\delta}^{2})^{2}\,w_{3}^{2}+(w_{1}^{2}-{\delta}^{2})^{2}(w_{2}^{2}-{\delta}^{2})\,w_{0}^{2}+(2w_{1}w_{2}+\bar{\nu})(w_{1}^{2}-{\delta}^{2})(w_{2}^{2}-{\delta}^{2})\,w_{3}w_{0}\\\
&+&\big{(}2w_{1}^{3}w_{2}^{2}+\bar{\nu}w_{1}^{2}w_{2}+\bar{c}w_{1}+\bar{a}-(2w_{1}w_{2}^{2}+\bar{\nu}w_{2}){\delta}^{2}-w_{1}{\delta}^{4}\big{)}(w_{2}^{2}-{\delta}^{2})\,w_{3}\\\
&+&\big{(}2w_{1}^{2}w_{2}^{3}+\bar{\nu}w_{1}w_{2}^{2}+\bar{c}w_{2}+\bar{a}-(2w_{1}^{2}w_{2}+\bar{\nu}w_{1}){\delta}^{2}-w_{2}{\delta}^{4}\big{)}(w_{1}^{2}-{\delta}^{2})\,w_{0}\\\
&+&w_{1}^{4}w_{2}^{4}+\bar{\nu}w_{1}^{3}w_{2}^{3}+\bar{c}w_{1}^{2}w_{2}^{2}+\bar{a}w_{1}w_{2}(w_{1}+w_{2})\\\
&-&\Big{(}\big{(}w_{1}^{2}w_{2}^{2}+\bar{\nu}w_{1}w_{2}\big{)}(w_{1}^{2}+w_{2}^{2})+\bar{a}(w_{1}+w_{2})\Big{)}{\delta}^{2}+(w_{1}^{2}w_{2}^{2}+\bar{\nu}w_{1}w_{2}-\bar{c}){\delta}^{4}-(w_{1}^{2}+w_{2}^{2}){\delta}^{6};\end{array}$
(5.3)

so the map (5.1) is Liouville integrable.

The Hamiltonian vector field for the first flow, generated by $H_{1}$, takes
the form

$\frac{\mathrm{d}w_{n}}{\mathrm{d}t}=(w_{n}^{2}-{\delta}^{2})(w_{n+1}-w_{n-1})$
(5.4)

for $n=1,2$, and once again, since the Poisson map (5.1) is compatible with
the flow $\\{\cdot,H_{1}\\}$, the equation (5.4) holds for all
$n\in{\mathbb{Z}}$, and thus the map and the flow together produce a sequence
of functions $\big{(}w_{n}(t)\big{)}_{n\in{\mathbb{Z}}}$ satisfying (5.4),
which (up to rescaling) is the general form of the modified Volterra lattice
equation. If we set ${\delta}\to 0$ in (5.4), then the equation (4.4) is
recovered, corresponding to the same limit that reproduces (4.1) as a
degenerate case of (5.1). However, the behaviour of the degenerate map (4.1)
is sufficiently different compared with (5.1) e.g. with respect to singularity
structure, that it is worth giving it a separate analysis as we have done
here.

Let us denote a solution of the Volterra lattice (3.6) by $\hat{w}_{n}$. Then
the Miura map from the modified Volterra lattice (5.4) takes the form

$\hat{w}_{n}=(w_{n+1}\mp{\delta})(w_{n}\pm{\delta}),$ (5.5)

(so there are effectively two maps, with an opposite choice of sign in each
factor on the right-hand side above). Moreover, this persists at the level of
the maps (3.2) and (5.1), in the following sense.

######  Theorem 5.1.

Let $w_{n}$ be a solution of (5.1) with parameters
$\bar{a},\bar{c},\bar{\nu}$, lying on the level set $H_{1}=\bar{h}_{1}$,
$H_{2}=\bar{h}_{2}$ of the first integrals (5.2) and (5.3). Then for either
choice of signs, $\hat{w}_{n}$ given by the Miura map (5.5) is a solution of
(3.2) with parameters

$\nu=\bar{\nu}+6{\delta}^{2},\qquad
b=\bar{c}+4\bar{\nu}{\delta}^{2}+7{\delta}^{4},\qquad
a=\bar{h}_{1}+\bar{c}{\delta}^{2}+\bar{\nu}{\delta}^{4}-{\delta}^{6}.$

Moreover, on either solution $\hat{w}_{n}$, the values $h_{1},h_{2}$ of the
first integrals (3.3) and (3.4) for the map (3.2) are given by

$h_{1}=\bar{h}_{2}+2{\delta}^{8},\qquad
h_{2}=-\bar{a}^{2}-\bar{\nu}\bar{h}_{2}-\bar{c}\bar{h}_{1}-2\bar{h}_{2}{\delta}^{2}+(\bar{h}_{1}-\bar{\nu}\bar{c}){\delta}^{4}-\bar{\nu}{\delta}^{8}-4{\delta}^{10}.$

  * Proof:

The first part of this result is verified by substituting the Miura formula
(5.5) directly into (3.2), using (5.1) to eliminate $w_{n+5}$ followed by
$w_{n+4}$, and then using (5.2) to eliminate $w_{n+3}$ on the level set
$H_{1}=\bar{h}_{1}$. After the initial substitution of the Miura map and
eliminating, all the final results are quadratic in ${\delta}$, so do not
depend on the choice of sign in (5.5). Similar calculations using the same
substitutions in the formulae (3.3) and (3.4), together with the expression
(5.3) on the level set $H_{2}=\bar{h}_{2}$, produce the expressions for
$h_{1},h_{2}$, which are the corresponding values of the first integrals for
(3.2). ∎

We can also malke use of a tau function substitution for (P.vi), which has the
more complicated structure

$w_{n}+{\delta}=\rho_{n}\,\frac{\sigma_{n+2}\tau_{n}}{\sigma_{n+1}\tau_{n+1}},$
(5.6)
$w_{n}-{\delta}=\frac{1}{\rho_{n+1}}\,\frac{\sigma_{n}\tau_{n+2}}{\sigma_{n+1}\tau_{n+1}},$
(5.7)

with

$\rho_{n+2}=\rho_{n}.$

This implies that

$\hat{w}_{n}^{(+)}=(w_{n}-{\delta})(w_{n+1}+{\delta})=\frac{\sigma_{n}\sigma_{n+3}}{\sigma_{n+1}\sigma_{n+2}},$
(5.8)
$\hat{w}_{n}^{(-)}=(w_{n}+{\delta})(w_{n+1}-{\delta})=\frac{\tau_{n}\tau_{n+3}}{\tau_{n+1}\tau_{n+2}}$
(5.9)

are both solutions of (3.2), and both sequences $(\sigma_{n})$ and
$(\tau_{n})$ satisfy the same Somos-9 relation. Thus the two different
formulae for the Miura map in (5.5) can be regarded as defining a Bäcklund
transformation for the discrete equation (3.2) with parameter ${\delta}$,
since given $\hat{w}_{n}^{(-)}$ and a solution $w_{n}$ of (5.1), a new
solution $\hat{w}_{n}^{(+)}$ of (3.2) is generated by taking

$\hat{w}_{n}^{(+)}=\hat{w}_{n}^{(-)}+2{\delta}(w_{n+1}-w_{n}).$

##  6 Conclusion

We have shown that the integrable maps (P.iv), (P.v) and (P.vi) from [7] are
closely related to one another, via Miura-type transformations, and they
provide genus two solutions of Volterra and modified Volterra lattices,
respectively. So far we do not have a complete understanding of what the
relations between these maps mean geometrically, particularly from the Poisson
and algebro-geometric points of view. However, since the construction of the
integrable maps ${\cal V}_{g}$ presented in [10] is valid for any $g\geq 1$,
this strongly suggests that (P.v) and (P.vi) should each be the $g=2$ members
of a family of maps defined for any $g$. In the elliptic case ($g=1$) we have
constructed elliptic solutions of the modified Volterra and Volterra lattices,
and showed how they are linked by the Miura transformation, essentially
recovering the solutions found in [24], which can be interpreted in terms of
integrable maps in the plane (QRT type). The complete description of these
results, together with the proposed extension to families of maps for all
$g\geq 1$, is planned for future work.

Acknowledgments: The research of ANWH was supported by Fellowship EP/M004333/1
from the Engineering & Physical Sciences Research Council, UK, extended by
EP/V520718/1 COVID 19 Grant Extension Allocation University of Kent, and the
grant IEC\R3\193024 from the Royal Society; he is also grateful to the School
of Mathematics and Statistics, University of New South Wales, for hosting him
during 2017-2019 as a Visiting Professorial Fellow with funding from the
Distinguished Researcher Visitor Scheme, and to Wolfgang Schief, who provided
additional support during his time in Sydney. ANWH also thanks DICATAM for
supporting his visit to Brescia in November 2022. FZ acknowledges the support
of Università di Brescia, GNFM-INdAM and INFN, Sezione di Milano-Bicocca, Gr.
IV - Mathematical Methods in NonLinear Physics (Milano, Italy).

## References

  * [1] M. Adler, P. van Moerbeke, and P. Vanhaecke.  Algebraic integrability, Painlevé geometry and Lie algebras, vol. 47 of Ergebnisse der Mathematik und ihrer Grenzgebiete. Springer-Verlag, Berlin, 2004.
  * [2] M. Bruschi, O. Ragnisco, P.M. Santini and T. Gui-Zhang, Integrable symplectic maps, Physica D 49 (1991) 273–294.
  * [3] X.-K. Chang, X.-B. Hu and G. Xin, Hankel Determinant Solutions to Several Discrete Integrable Systems and the Laurent Property, SIAM J. Discrete Math. 29 (2015) 667–682.
  * [4] X.-M. Chen, X.-B. Hu and F. Müller-Hoissen, Non-isospectral extension of the Volterra lattice hierarchy, and Hankel determinants, Nonlinearity 31 (2018) 4393–4422.
  * [5] C. Cresswell and N. Joshi, The discrete first, second and thirty-fourth Painlevé hierarchies, J. Phys. A: Math. Gen. 32 (1999) 655–69.
  * [6] G. Gubbiotti, N. Joshi, D.T. Tran and C.-M. Viallet, Complexity and integrability in 4D bi-rational maps with two invariants. In: Nijhoff F., Shi Y., Zhang D. (eds), Asymptotic, Algebraic and Geometric Aspects of Integrable Systems. Springer Proceedings in Mathematics & Statistics, vol. 338, Springer, Cham, pp.17–36, 2020.
  * [7] G. Gubbiotti, N. Joshi, D.T. Tran and C.-M. Viallet, Bi-rational maps in four dimensions with two invariants, J. Phys. A: Math. Theor. 53 (2020) 115201.
  * [8] G. Gubbiotti, Lagrangians and integrability for additive fourth-order difference equations, Eur. Phys. J. Plus 135 (2020) 853.
  * [9] A.N.W. Hone, Continued fractions and Hankel determinants from hyperelliptic curves, Comm. Pure Appl. Math. 73 (2020) 2310–2347.
  * [10] A.N.W. Hone, J.A.G. Roberts and P. Vanhaecke, A family of integrable maps associated with the Volterra lattice, arXiv:2309.02336
  * [11] N. Joshi and C.-M. Viallet, Rational maps with invariant surfaces, J. Int. Syst. 3 (2018) xyy017.
  * [12] M. Kac and P. van Moerbeke, On an Explicitly Soluble System of Nonlinear Differential Equations Related to Certain Toda Lattices, Adv. Math. 16 (1975) 160–169.
  * [13] D. Levi and R. Benguria, Bäcklund transformations and nonlinear differential difference equations, Proc. Natl. Acad.Sci. USA 77 (1980) 5025–5027.
  * [14] D. Levi and R. Yamilov, Conditions for the existence of higher symmetries of evolutionary equations on the lattice, J. Math. Phys. 39 (1997) 6648–6674.
  * [15] D. Levi, Bäcklund transformations, in Encyclopedia of Mathematical Physics, pp. 241–247, Academic Press, 2006.
  * [16] S. Maeda, Completely integrable symplectic mapping, Proc. Japan Acad. Ser. A Math. Sci. 63 (1987) 198–200.
  * [17] G.R.W. Quispel, H.W. Capel and J.A.G. Roberts, Duality for discrete integrable systems, J. Phys. A: Math. Gen. 38 (2005) 3965–3980.
  * [18] J. Shohat, On Stieltjes Continued Fractions, Amer. J. Math. 54 (1932) 79–84.
  * [19] T.-J. Stieltjes, Recherches sur les fractions continues, Ann. Fac. Sci. Toulouse Math. 8 (1894) J1-J122.
  * [20] Y.B. Suris, The Problem of Integrable Discretization: Hamiltonian Approach, Birkhäuser, 2003.
  * [21] P. Vanhaecke, Integrable systems in the realm of algebraic geometry, 2nd edition, Springer, 2001.
  * [22] A.P. Veselov, Integrable Maps. Russ. Math. Surv. 46 (1991) 1–51.
  * [23] D. Wójcik and J. Cieśliński (eds.), Nonlinearity & Geometry. Proceedings of the First Nonorthodox School on Nonlinearity and Geometry, Polish Scientific Publishers PWN, Warszawa, 1998.
  * [24] Z. Yan, Discrete exact solutions of modified Volterra and Volterra lattice equations via the new discrete sine-Gordon expansion algorithm, Nonlinear Analysis 64 (2006) 1798–1811.
  * [25] J.-P. Yu, W.-X. Ma, Y.-L. Sun and C.M. Khalique, N-fold Darboux transformation and conservation laws of the modified Volterra lattiice, Mod. Phys. Lett. B 32 (2018) 1850409.