# Nonequilibrium study of the $J_{1}-J_{2}$ Ising model with random $J_{2}$
couplings in the square lattice

Octavio D. Rodriguez Salmon  Departamento de Física, Universidade Federal do
Amazonas, 3000, Japiim, 69077-000, Manaus-AM, Brazil<EMAIL_ADDRESS>Minos A. Neto  Departamento de Física, Universidade Federal do Amazonas, 3000,
Japiim, 69077-000, Manaus-AM, Brazil  Thiago Lobo  Instituto Federal de
Educação, Ciência e Tecnologia do Rio de Janeiro - IFRJ
Rua Dr. José Augusto Pereira dos Santos, s/n (C.I.E.P. 436 Neusa Brizola),
24425-004 São Gonçalo, Rio de Janeiro, Brasil.  Francisco Dinóla Neto  Centro
Universitário do Norte - UNINORTE
Rua Huascar de Figueiredo, 290, Centro
69020-220, Manaus, AM, Brazil.

###### Abstract

We studied the critical behavior of the $J_{1}-J_{2}$ spin-1/2 Ising model in
the square lattice by considering $J_{1}$ fixed and $J_{2}$ as random
interactions following discrete and continuous probability distribution
functions. The configuration of $J_{2}$ in the lattice evolves in time through
a competing kinetics using Monte Carlo simulations leading to a steady state
without reaching the free-energy minimization. However, the resulting non-
equilibrium phase diagrams are, in general, qualitatively similar to those
obtained with quenched randomness at equilibrium in past works. Accordingly,
through this dynamics the essential critical behavior at finite temperatures
can be grasped for this model. The advantage is that simulations spend less
computational resources, since the system does not need to be replicated or
equilibrated with Parallel Tempering. A special attention was given for the
value of the amplitude of the correlation length at the critical point of the
superantiferromagnetic-paramagnetic transition.

##  I Introduction

Nonequilibrium classical spin Hamiltonians with random parameters are still
attracting the studies of Statistical Mechanics mijat . Problems like spin
glasses were generally studied by spin models with quenched randomness, where
the couplings between pairs of spins are fixed in time, but their values are
randomly distributed throughout the lattice. So this cause frustration, which
is a fundamental ingredient for this glassy state. In this context the model
is studied at equilibrium. Nevertheless, this neglects the possibility that
frustration could avoid the system to reach an equilibrium state and so the
macroscopic behavior would be determined by the dynamics.

In order to simulate this situation, a competing kinetic process was taken
into account to induce the presence of nonequilibrium steady states
garrido_1991 . On the other hand, quenched randomness disregard diffusion of
magnetic ions lopez_1992 . When diffusion occurs, the distance between pairs
of spins varies with time. A way of modeling this phenomena in a lattice,
consists in considering the exchange couplings $J_{ij}$ to change in time and
space gonzalez_1994 . Past studies reported that the diffusion of disorder
affect the macroscopic behavior of the system gonzalez_1994 .

A kind of conflicting dynamics was implemented to study the random-field Ising
model (RFIM) nuno_2009 ; nuno_2010  inspired in the work of Gonzalez et. al.
gonzalez_1994 . Interestingly, it has been demonstrated paula_2014  that
quenched randomness destroys any first-order phase transition in low
dimensions. Thus, the RFIM does not present a tricritical point in the
ferromagnetic-paramagnetic frontier bellow the upper critical dimension
($d_{c}=6$) hartman_2011 . However, in the context of this competing kinetic,
the phase diagram of the RFIM exhibits a tricritical point when the field
obeys a double-Gaussian probability distribution in the square lattice
nuno_2009 , and for a bimodal distribution in the cubic lattice nuno_2010 .
Accordingly, these facts motivate the comparison between results with quenched
randomness and this conficting dynamics.

In this paper we studied the $J_{1}-J_{2}$ Ising Model with in the particular
case where the next-nearest-neighbor interactions $J_{2}$ are random, and
evolve in time by a competing kinetics which is detailed in the next section
nuno_2009 ; nuno_2010 . This model with fixed $J_{1}$ and $J_{2}$
antiferromagnetic interactions has shown controversial results in its critical
behavior nightingale_1977 ; swendsen_1979 ; otimaa_1981 ; binder_1980 ;
landau_1980 ; landau_1985 ; landau_2000 . Its phase diagram is studied in the
plane where the horizontal axis is either $r=|J_{1}/J_{2}|$ or
$r=|J_{2}/J_{1}|$ , and the vertical axis is the temperature $T$. For
$J_{2}<J_{1}$ the low temperature order is antiferromagnetic AF, and for
$J_{2}<J_{1}$ the colinear or superantiferromagnetic SAF order is present. It
was believed that for $J_{1}/2<J_{2}\lesssim J_{1}$ the critical frontier was
continuous with weak universality, i.e. exponents like $\beta$ and $\gamma$
break the universality, but their normalized exponents $\beta/\nu$ and
$\gamma/\nu$ remain the same. However, a later work showed a double-peaked
structure of the energy histogram suggesting weak first-order transition
kalz_2011 .

When $J_{2}$ takes the value $\lambda_{1}J_{1}$ with probability $p$ and
$-\lambda_{2}J_{1}$ with probability $1-p$, a competition emerges between SAF
and AF phases if $J_{1}$ is antiferromagnetic, and between SAF and F if
$J_{1}$ is ferromagnetic (with $\lambda_{1}>0$ and $\lambda_{2}\geq 0$). These
two cases were studied considering quenched randomness. In the first one the
authors simulated the model with bond dilution in $J_{2}$ ($\lambda_{1}=1$ and
$\lambda_{2}=0$) yining_2018 . They used Parallel Tempering so as to
equilibrate the system. This Monte Carlo technique is commonly used for
equilibrate models with quenched randomness at low temperatures. So it must be
executed in parallel for simulating each replica of the system at different
temperatures for each bond realization swendsen_1986 ; david_2005 .

The authors yining_2018  obtained a phase diagram in the plane $T$ versus the
degree of dilution $x=p$, where the frontiers SAF-P and AF-P are separated by
a gap between $0.3<x<0.8$ at $T=0$. They found strong evidences of a spin-
glass phase in this gap only at zero temperature, regarding that there should
not be spin-glass phase in two dimensions for Ising-type spin-glass models
(see references nobre_2001  and fernandez_2016  for this controversy). On the
other hand, the second case with $\lambda_{1}=\lambda_{2}=1$ was treated in
the effective-field approximation obtaining a phase diagram with similar
topological characteristics octavio_2009 , though the nature of the gap at
$T=0$ remains an open question. So, in this work we mainly simulate these two
cases in which the random bonds $J_{2}$ change in time by a conflicting
dynamics nuno_2009 ; nuno_2010 .

##  II The Model and Simulation

The system is modeled by the following Hamiltonian:

$\displaystyle{\cal
H}(\boldsymbol{S,J_{2}})=-J_{1}\sum_{(nn)}S_{i}S_{j}-\sum_{(nnn)}J_{2}^{(i,j)}S_{i}S_{j},$
(1)

where $S_{i}=\pm 1$, for $i=1,2,3,...,N$, with $N=L\times L$ is the number of
sites of a square lattice of size $L$. $\mathbf{S}=\\{S_{1},S_{2},\dots
S_{N}\\}$ is the set of spin variables, and
$\mathbf{J_{2}}=\\{J_{2}^{(i,j)}\\}$ is the set next-nearest-neighbor coupling
values.

In this work the following probability distributions functions (PDF) for
$J_{2}$ are used:

$\displaystyle{\cal
P}(J_{2}^{(i,j)})=p\delta(J_{2}^{(i,j)}-J_{1})+(1-p)\delta(J_{2}^{(i,j)}+J_{1}),$
(2)

where $0\leq p\leq 1$ and $J_{1}>0$ (ferromagnetic nearest-neighbor
couplings). So by this PDF each coupling $J_{2}^{(i,j)}$ takes the value
$J_{1}$ with probability $p$, and the value $-J_{1}$ with probability $1-p$.
This is the case studied with quenched randomness in reference octavio_2009 .
Then

$\displaystyle{\cal P}(J_{2}^{(i,j)})=p{\cal
G}(J_{2}^{(i,j)}-J_{1})+(1-p){\cal G}(J_{2}^{(i,j)}+J_{1}),$ (3)

where $J_{1}>0$ and the function is a Gaussian one:

$\displaystyle{\cal
G}(u)=\frac{1}{\sqrt{2\pi}\sigma}\exp(-\frac{u^{2}}{2\sigma^{2}}),$ (4)

where $\sigma$ is the standard deviation. Then

$\displaystyle{\cal
P}(J_{2}^{(i,j)})=p\delta(J_{2}^{(i,j)})+(1-p)\delta(J_{2}^{(i,j)}-J_{1}),$
(5)

where $J_{1}<0$ (antiferromagnetic nearest-neighbor couplings). In this PDF
each coupling $J_{2}^{(i,j)}$ is zero with probability $p$ (bond dilution),
and takes the value $J_{1}$ with probability $1-p$. This is the PDF studied
with quenched randomness in reference yining_2018 . The last distribution used
in this work is this double Gaussian PDF :

$\displaystyle{\cal P}(J_{2}^{(i,j)})=p{\cal G}(J_{2}^{(i,j)})+(1-p){\cal
G}(J_{2}^{(i,j)}-J_{1}),$ (6)

where $J_{1}<0$. Note that for $\sigma\to 0$, we recover the bimodal
distribution given in Eq.(5). The system is in contact with a thermal
reservoir, which provides the temperature $T$ of the system inducing
stochastic changes given by the following master equation:

$\displaystyle\frac{\partial P(t,\boldsymbol{S^{b}})}{\partial
t}=\sum_{\boldsymbol{S^{a}}}(W(\boldsymbol{S^{a}}\rightarrow\boldsymbol{S^{b}})P(t,\boldsymbol{S^{a}})-W(\boldsymbol{S^{b}}\rightarrow\boldsymbol{S^{a}})P(t,\boldsymbol{S^{b}})),$
(7)

being $W(\boldsymbol{S^{a}}\rightarrow\boldsymbol{S^{b}})$ the probability per
unit time for a transition from the spin configuration $\boldsymbol{S^{a}}$ to
$\boldsymbol{S^{a}}$, and $P(t,\boldsymbol{S^{b}})$ the probability that the
system is in the spin configuration $\boldsymbol{S^{b}}$ at time $t$. A
stationary solution of the master equation is given by the detailed balance
condition
$W(\boldsymbol{S^{b}}\rightarrow\boldsymbol{S^{a}})=\exp(-\beta\Delta{\cal
H})$, where $\Delta{\cal H}={\cal H}(\boldsymbol{S^{a},J_{2}})-{\cal
H}(\boldsymbol{S^{b},J_{2}})$, and $\beta=1/(k_{B}T)$. So, a steady state is
reached by the Metropolis Algorithm
$W(\boldsymbol{S^{b}}\rightarrow\boldsymbol{S^{a}})=\min\\{1,\exp(-\beta\Delta{\cal
H})\\}$ gonzalez_1994 ; nuno_2009 ; nuno_2010 .

To implement it in a simulation a Monte Carlo sweep (MCS) consists in
generating a new configuration of random $J_{2}$ bonds according to a given
PDF; then, all lattice sites are visited, and each spin ﬂips according to the
Metropolis’ criterium. In this way the random bonds $J_{2}$ vary with time by
changing their values after each MCS. This is the essence of this conflict
dynamics where there is a competing kinetic between spins and random bonds.
Thus the system reaches a nonequilibrium steady state. This context differs
from those of equilibrium quenched and annealed random models.

For a given lattice size $L$ of the system, we note that the equilibration is
guaranteed after $10^{5}$ MCS for a given temperature. To overcome the
critical slowing down we accumulated values each $L$ MCS. In some cases we run
up to $3.6\times 10^{6}$ MCSs to reduce fluctuations in the thermal mean about
the critical temperature. In this way the order parameter, the susceptibility,
the Binder Cumulant and the correlation length were averaged for a given
temperature. In this work we used three order parameters according to the
phase formed by the parameters established in a particular simulation. The
first one is the magnetization per spin, which measures the ferromagnetic
order F:

$\displaystyle m_{F}=|\frac{1}{N}\sum_{i=1}^{N}S_{i}|.$ (8)

The antiferromagnetic order parameter meassuring the Neel order or the AF
order is as follows :

$\displaystyle m_{AF}=|\frac{1}{N}\sum_{i=1}^{N}(-1)^{x_{i}+y_{i}}S_{i}|.$ (9)

The order parameter for the colinear order (or the superantiferromagnetic
order) is given by :

$\displaystyle m_{SAF}$ $\displaystyle=\sqrt{m_{x}^{2}+m_{y}^{2}}$ (10)
$\displaystyle m_{x}$
$\displaystyle=\frac{1}{N}\sum_{i=1}^{N}(-1)^{x_{i}}S_{i}$ $\displaystyle
m_{y}$ $\displaystyle=\frac{1}{N}\sum_{i=1}^{N}(-1)^{y_{i}}S_{i}.$

The susceptibility corresponding to each of these order parameters is given
by:

$\displaystyle\chi=\frac{\langle\phi^{2}\rangle-{\langle\phi\rangle}^{2}}{k_{B}T},$
(11)

where $\phi$ is any of the above order parameters and $\langle\dots\rangle$ is
the thermal mean, which is the average calculated from all the meassurements
taken after the steady state is reached, for a given temperature. On the other
hand, the Binder cumulant is useful for estimating the critical points, and is
calculated as follows:

$\displaystyle
g=1-\frac{\langle\phi^{4}\rangle}{3{\langle\phi^{2}\rangle}^{2}}$ (12)

As is well known, the crossing point of $g_{L}-T$ curves for different sizes
$L$ locates the critical point for a second-order phase transition. For a
first-order transition the cummulant presents a particular shape with a low
peak at which the transition temperature can be localized. In this work we did
not find evidences of a first-order behavior. Furthermore, we also used the
second moment correlation length $\xi$ to locate the critical point
katzgraber2008 ; campbell_2019 , which is given by:

$\displaystyle\xi$
$\displaystyle=\frac{1}{2\sin(k_{min}/2)}\sqrt{\frac{\widetilde{G}(\boldsymbol{0})}{\widetilde{G}(\boldsymbol{k_{min}})}-1}$
(13) $\displaystyle\widetilde{G}(\boldsymbol{k})$
$\displaystyle=\frac{1}{N}\sum_{\boldsymbol{r}}G(\boldsymbol{r})\exp(i\boldsymbol{k.r}),$

where $G(\boldsymbol{r})$ the correlation function, $\boldsymbol{r}=(x,y)$ and
$\boldsymbol{k_{min}}=(2\pi/L,0)$. In order to associate the correlation with
the leading magnetic orderings studied in this work, this formula is applied
in the following form:

$\displaystyle\xi$
$\displaystyle=\frac{1}{2\sin(\pi/L)}\sqrt{\frac{\langle\phi^{2}\rangle}{\langle
F\rangle}-1},$ (14) $\displaystyle F$
$\displaystyle=\frac{1}{2N}(F_{1}+F_{2}),$

where $F_{1}$ and $F_{2}$ are computed according to the corresponding order
parameter, as follows

$\displaystyle F_{1}$
$\displaystyle=\left(\sum_{x=1}^{L}\sum_{y=1}^{L}a_{x,y}S_{x,y}\cos(2\pi
x/L)\right)^{2}+\left(\sum_{x=1}^{L}\sum_{y=1}^{L}a_{x,y}S_{x,y}\sin(2\pi
x/L)\right)^{2}$ (15) $\displaystyle F_{2}$
$\displaystyle=\left(\sum_{x=1}^{L}\sum_{y=1}^{L}b_{x,y}S_{x,y}\cos(2\pi
y/L)\right)^{2}+\left(\sum_{x=1}^{L}\sum_{y=1}^{L}b_{x,y}S_{x,y}\sin(2\pi
y/L)\right)^{2}.$

If $\phi=m_{F}$, $\xi$ is denoted by $\xi_{F}$, and $a_{x,y}=b_{x,y}=1$. If
$\phi=m_{AF}$, $\xi=\xi_{AF}$ and $a_{x,y}=b_{x,y}=(-1)^{x+y}$. When
$\phi=m_{SAF}$, $\xi=\xi_{SAF}$, $a_{x,y}=(-1)^{x}$ and $b_{x,y}=(-1)^{y}$.

The formula of Eq.(14) is valid for temperatures greater or equal to the
critical temperature, where $\langle\phi\rangle=0$. In the same way the
crossing point of curves with different sizes of the Binder Cummulant locates
the critical point, different curves of $\xi/L$ are intercepted at the
transition point of a second-order phase transition (unless scaling
corrections were needed). For instance, the correlation length defined in
Eq.(14) depends on the lattice size $L$ at $T=T_{c}$, as follows:

$\xi=L(x^{*}+AL^{-\Delta}+\dots).$ (16)

So, at $T=T_{c}$, the limit $\lim_{L\to\infty}\xi(L)/L$ tends to a constant,
which is $x^{*}\approx 0.905$, for the 2D Ising model salas2000 . On the other
hand, at higher dimensions, above the upper critical dimension, the scaling
behavior is different from that of Eq.(16) flores2015 ; young2005 .

##  III Results

In what follows, the value of $J_{2}$ will be expressed in units of $J_{1}$,
and the critical temperature in $J_{1}/k_{B}$ units. Thus, these values are
set as 1.

###  III.1 Case with ferromagnetic nearest-neighbor interactions ($J_{1}>0$)

For this case we considered firstly the $J_{2}$ bonds obeying the PDF given in
Eq (2). So we have a spin lattice with ferromagnetic nearest-neighbor
interactions with random bonds being ferromagnetic with probability $p$ and
superantiferromagnetic with probability $1-p$. An exploration of the order
parameters $m_{F}$ and $m_{SAF}$ given in equations (8) and (10) was done for
$0\leq p<1$ so as to obtain the frontiers dividing phases SAF and F with the
paramagnetic one (P).

For $p=1$, we have the trivial case of the Ising model in the square lattice
with nearest- and next-nearest-neighbor interactions having a second-order
phase transition at $T_{c}=5.262\pm 0.001$ (in agreement with reference
dalton_1969 ). For $p=0$, the SAF phase is present at low temperatures, and
the SAF-P transition is found at $T=2.083\pm 0.002$ (the same value was
obtained in malakis_2006 ). Figures 1 and 2 show how these values were
determined by locating the crossing point of respective Binder Cummulant and
correlation length. It is important to point out that the present curves of
the second moment correlation length are only valid for $T\geq T_{c}$, but for
the sake of visualization they are plotted for $T<T_{c}$ too. Of course, the
criticality for $p=0$ and $p=1$ is reached at equilibrium, since no competing
kinetics takes place.

The role of $p$, for $0<p<1$, is to create a competition between SAF and F,
which is stronger for its intermediate values (see Eq.(2)). When $p$ is close
to 1, the ferromagnetic order is present as Fig.3 shows. There it is observed
the Binder Cummulant $g_{F}$ and the second moment correlation length divided
by the lattice size $\xi_{F}/L$ as functions of the temperature, corresponding
to the ferromagnetic order parameter $m_{F}$, for $p=0.6$. Four curves, for
sizes $L=30,60,90,120$, colapse at a point in Fig.3(a) and (b) locating the
critical temperature. Both figures agree with the value $T_{c}=1.410\pm
0.005$. This is a clear signal of second-order phase transition. However, the
scale invariance at the critical point seems to be better visualized using the
correlation length.

In Fig.4 is shown the correlation length $\xi_{SAF}$ and $\xi_{F}$ scanned
through the probability $p$, for a very low temperature $T=0.05$. In Fig.4a
the crossing point is located at $p_{c}=0.0783\pm 0.0004$, whereas in Fig.4b
we have $p_{c}=0.571\pm 0.006$. Consequently, there is a paramagnetic gap
between the SAF-P and F-P second-order frontiers. In Fig.5 the resulting phase
diagram is exhibited, where all points were obtained using the correlation
length and the Binder Cummulant. Note the paramagnetic gap for $p_{c1}\leq
p\leq p_{c2}$, where $p_{c1}\approx 0.078$ and $p_{c2}\approx 0.57$.

If $J_{2}$ obeys the PDF given in Eq.(3) (a double Gaussian), any finite value
of the standard deviation $\sigma$ destroys the SAF phase. In what follows,
$\sigma$ will be expressed in units of $J_{1}$. Note that the PDF given in
Eq.(2) is the limit of the double Gaussian given in Eq.(2), when $\sigma\to
0$. Interestingly, the effect of a very small $\sigma$ is to shift the extreme
critical points of the phase diagram from $p_{c}\approx 0.57$ to $p_{c}\approx
1/3$ and from $T_{c}(p=1)\approx 5.26$ to $T_{c}(p=1)\approx 3.57$. We
confirmed that this shift exists even for a very small value of $\sigma$, such
as $\sigma=0.001$.

In Fig.5 the phase diagram is presented for $\sigma=0.1$. There, the extreme
points of the frontier are $p_{c}=0.3310\pm 0.0005$, for $T=0.005$, and for
$p=1$, $T_{c}=3.5730\pm 0.0004$. To understand the effect of the standard
deviation on $p_{c}$ we show in Fig7a how $p_{c}$ evolves with $\sigma$. It
remains constant for $0<\sigma<0.3$, because the gaussians of the PDF do not
overlap, as Fig7b shows. Then, it increases from $\sigma\approx 0.3$, since
the gaussians of the PDF increases their overlap. For $\sigma=2.5$,
$p_{c}=0.987$(see Fig7b). Accordingly, for $\sigma\gtrsim 2.5$, $p_{c}$
reaches its maximum value, which is one. So, the ferromagnetic area of the
phase diagram is totally reduced, since the randomness caused by $\sigma$ is
so strong, such that any amount of the second Gaussian in Eq.(3) (any amount
of ($1-p$)), would destroy the ferromagnetic phase (see Fig7). On the other
hand, the effect of $\sigma$ on $T_{c}(p=1)$ is presented in Fig.8.

In Fig.8a it can be observed the discontinuous fall in the critical
temperature $T_{c}(p=1)$ when $\sigma$ assumes infinitesimal values, remaining
constant for $0<\sigma<0.5$. Then it increases slightly before reaching its
maximum for $\sigma\approx 1$. Then $T_{c}(p=1)$ decreases. The extension of
the Gaussian distribution for these three different stages is exhibited in
Fig.8b, for three representative values of $\sigma$. Note that the critical
temperature is unaltered when the probability for negative values of $J_{2}$
is zero. However, if the left tail of the Gaussian includes the interval
$-2\lesssim J_{2}<0$, $T_{c}(p=1)$ increases. If this interval of negative
$J_{2}$ values is extended, the effect of $\sigma$ is to decrease the critical
temperature. In Fig.8a the critical temperature is plotted up to $\sigma=2.5$,
when $p_{c}$ assumes the value 1.

###  III.2 Case where $J_{1}<0$

Now we treat the case $J_{1}<0$, where the nearest-neighbor interactions are
antiferromagnetic. If the PDF given in Eq.(5) is obeyed by the $J_{2}$ bonds,
these are antiferromagnetic with probability $1-p$, and zero with probability
$p$. So, $p$ measures the degree of dilution of $J_{2}$ bonds. This particular
case was studied in equilibrium with quenched randomness in reference
yining_2018 , where $x$ stood for our probability of dilution $p$, as
mentioned before in the Introduction.

The phase diagram we obtained is shown in Fig.9, where the SAF-P and AF-P
frontiers are separated by paramagnetic gap for $p_{c1}\leq p\leq p_{c2}$,
where $p_{c1}\approx 0.309$ and $p_{c2}\approx 0.69$. The critical
temperatures at the extremes of the phase diagram are $T_{c}(p=0)\approx 2.08$
(in agreement with malakis_2006 ), and $T_{c}(p=1)=2.269...$ (the exact 2D
Ising one).

Phase diagram in Fig.9 can be compared with that exhibited in Fig.5a of
reference yining_2018 . There, the Neel phase is our AF phase and the strip
antiferromagnetic one is our SAF phase. The topology is quite similar.
Nevertheless, we did not find a spin-glass phase in the gap area close to
$T=0$, as in reference yining_2018 . At least, not with the spin-glass order
parameter averaging the correlation between the spins of two copies of the
system (see Eq.(6) in yining_2018 ).

Furthermore, the frontiers we obtained exhibit humps, which do not appear with
quenched randomness in that paper. In the present case, the SAF-P frontier
reaches its maximum at $(p,T_{c})\approx(0.075,2.196)$, and at
$(p,T_{c})\approx(0.925,2.364)$, for the AF-P frontier, as shown in Fig.10.
These humps may be caused by the particular type of competing kinetic used in
this work, which causes an increase of the critical temperature for values of
$p$ close to zero and $1$.

Finally, we applied the PDF given in Eq.(6) for the $J_{2}$ bonds, with
$J_{1}<0$. In this case $J_{2}$ is distributed by two gaussians centered at
$J_{2}=0$ and $J_{2}=J_{1}$, with probability $p$ and $1-p$, respectively.
When $\sigma\to 0$ the bimodal given in Eq.(5) is recovered. In Fig11 is
presented the phase diagram for $\sigma=0.1$. Again we observe that the SAF is
not present. This happens even when $\sigma$ is infinitesimal. Thus, only the
AF survives in an area enclosed by a second-order frontier, where
$p_{c}\approx 0.692$ and $T_{c}(p=1)$ is the 2D Ising critical temperature.
The evolution of $p_{c}$ with $\sigma$ is shown in Fig.12a, where we may see
how $p_{c}$ falls from $p_{c}\approx 0.692$ down to $p_{c}\approx 0.383$, when
$\sigma$ assumes an infinitesimal value. Then, $\sigma$ does not affect the
value of $p_{c}$ until $\sigma\gtrsim 0.35$.

The AF area disappears when $p_{c}=1$, this happens for $\sigma\gtrsim 1.3$.
In Fig.12b is shown that $p_{c}$ is not affected by $\sigma$ when the tails of
the two gaussians die for $-1<J_{2}<1$. For $\sigma\lesssim 0.35$, the tails
of the two gaussians start to include $|J_{2}|>1$, so $p_{c}$ begins to
increase. Similarly, in Fig.13 is shown the evolution of $T_{c}(p=1)$ with
$\sigma$, until $p_{c}$ reaches the value 1. It confirms that the frontier
exhibited in Fig.11 is unaltered when $\sigma<0.35$.

###  III.3 The amplitude $x^{*}_{SAF}$

The amplitude $x^{*}_{SAF}$ in Eq.(16) has attracted our attention. We may
observe that the crossing point of $\xi_{SAF}/L$ curves gives different values
of $x^{*}_{SAF}$ when $J_{2}=-J_{1}<0$ and $J_{2}=J_{1}<0$. This is exhibited
in Fig.14, where the point of intersection provides $T_{c}(p=0)$ in Figures 5
and 9. In Fig.14a the uncertainty of the crossing point also gives
$x^{*}_{SAF}=1.210\pm 0.005$, and in Fig.14b $x^{*}_{SAF}=1.483\pm 0.008$.
However, the critical temperatures are just the same $T_{c}(p=0)\approx
2.082$, whithin the error.

So, though both critical points are just at the same SAF-P critical
temperature, the way $J_{1}$ and $J_{2}$ bonds are configured to form the SAF
phase gives a slight different value of $x^{*}_{SAF}$. To reinforce this idea,
we analysed the effect of the bimodal distributions given in Eqs.(2) and (5)
on this superanfiterromagnetic amplitude. For the case $J_{1}>0$, the SAF
phase is destroyed for a small value of $p$ (see Fig. 5), so no relevant
change in $x^{*}_{SAF}$ is observed. However, when $J_{1}<0$, the degree of
dilution of $J_{2}$ changes this amplitude, as suggested by an ad-hoc finite-
size criterium.

To show it, we have determined the interval of uncertainty of $x^{*}_{SAF}$ in
Fig.15, for $p=0$, using Eq.(16). We used small sizes of $L$, such as
$L=8,10,12,14,16,18,20,24,28$, for which the second term in Eq.(16) is
relevant. The limits of $x^{*}_{SAF}$ were determined by extrapolating to
$L\to\infty$ the line of points of $\xi_{SAF}/L$ versus $L^{-\Delta}$, each
one corresponding to the limits of the critical temperature found in the
crossing point of Fig.14b.

The values of $\Delta$ were scanned in order to get the adequate ones to
adjust the points in a straight line. So, the two extrapolated lines in Fig.15
intersect to the vertical axis at points which constitute the limits of the
uncertainty of $x^{*}_{SAF}$. Accordingly, we found $1.440<x^{*}_{SAF}<1.497$.
A similar extrapolation was performed for the transition point at $T=0.4$ of
the left frontier in Fig.9. There we have $0.30882<p_{c}<0.30898$, according
to the uncertainty of the crossing point shown in Fig. 16a.

The resulting interval of the amplitude is shown in Fig.16b. Thus,
$1.53<x^{*}_{SAF}<1.55$, which does not intersect the interval found for
$p=0$. Therefore, this finite-size analysis suggest that $x^{*}_{SAF}$ also
changes with the dilution of $J_{2}$. Of course, more rigorous numerical
treatments need to be done to confirm this affirmation.

##  IV Conclusions

We have studied the steady states of a nonequilibrium $J_{1}-J_{2}$ Ising
model with random $J_{2}$ bonds obeying bimodal and double Gaussian
probability distribution functions, for ferromagnetic and antiferromagnetic
$J_{1}$ couplings, in a competing kinetics.

All the phase diagrams presented second-order frontiers dividing the
paramagnetic phase (P) with some of the following orders : the F, AF or SAF
order. Other authors found similar topologies for these phase diagrams
obtained at equilibrium using quenched randomness with the bimodal
distributions considered in this work. We observed that the SAF phase can not
exist when $J_{2}$ obeys a continuous probability distribution.

Also, no universality breaking was observed for the critical exponents of the
corresponding order parameters and their susceptibilities. However, we found
that the amplitude $x^{*}_{SAF}$ of the correlation length at the SAF-P
transition is not only different to that of 2D pure Ising model, but may
change if the degree of dilution of $J_{2}$ bonds is relevant. This was shown
by calculating the uncertainties of $x^{*}_{SAF}$. As far as we know, no other
work has calculated this particular amplitude. It would be interesting to
perform the same calculations with quenched randomness in order to compare the
obtained $x^{*}_{SAF}$ values. This will be done in a future work.

ACKNOWLEDGEMENTS

This work was partially supported by CNPq - Grants: 408787/2018-0 and
306569/2018-3 (Brazilian Research Agency). We also thank Ivan Scivetti for
fruitful discussion.

## References

  * (1) Svetislav Mijatović, Dragutin Jovković, and Djordje Spasojević, Phys. Rev. E 103, 032147 (2021).
  * (2) P. L. Garrido and J. Marro, Europhys. Lett. 15, 375 (1991).
  * (3) A. I. López-Lacomba and J. Marro, Phys. Rev. B 46, 8244 (1992).
  * (4) J. M. Gonzalez-Miranda, A. Labarta, M. Puma, Julio F. Fernández, P. L Garrido and J. Marro, Phys. Rev. E 49 2041 (1994).
  * (5) N. Crokidakis, J. Stat. Mech. P02058 (2009).
  * (6) N. Crokidakis, Phys. Rev. E 81, 041138 (2010).
  * (7) P. Villa Martín, Juan A. Bonachela, and Miguel A. Muñoz. Phys. Rev. E 89, 012145 (2014).
  * (8) Björn Ahrens and Alexander K. Hartmann, Phys. Rev. B 83, 014205 (2011).
  * (9) M. P. Nightingale, Phys. Lett. A 59, 486 (1977).
  * (10) R. H. Swendsen and S. Krinsky, Phys. Rev. Lett. 43, 177 (1979).
  * (11) J. Oitmaa, J. Phys. A: Math. Gen. 14, 1159 (1981).
  * (12) K. Binder and D. P. Landau, Phys. Rev. B 21, 1941 (1980).
  * (13) D. P. Landau, Phys. Rev. B 21, 1285 (1980).
  * (14) D. P. Landau and K. Binder, Phys. Rev. B 31, 5946 (1985).
  * (15) D. P. Landau and K. Binder, Monte Carlo Simulations in Statistical Physics, 1st ed. (Cambridge University Press, Cambridge, 2000).
  * (16) A. Kalz, A. Honecker, and M. Moliner Phys. Rev. B 84, 174407, (2011).
  * (17) Yining Xu and Dao-Xin Yao Phys. Rev. B 97, 224419 (2018).
  * (18) R. H. Swendsen and J. S. Wang, Phys. Rev. Letters 57, 2607 (1986).
  * (19) David J. Earl and Michael W. Deem Phys. Chem. Chem. Phys. 7, 3910 (2005).
  * (20) F. D. Nobre, Phys. Rev. E 64, 046108 (2001).
  * (21) L. A. Fernandez, E. Marinari, V. Martin-Mayor, G. Parisi, and J. J. Ruiz-Lorenzo Phys. Rev. B 94, 024402 (2016).
  * (22) Octavio. D. R. Salmon, J. Ricardo de Sousa and F.D. Nobre, Phys. Lett. A, 373, 2525 (2009).
  * (23) Helmut G. Katzgraber, I. A. Campbell, and A. K. Hartmann, Phys. Rev. B 78, 184409 (2008).
  * (24) Ian A. Campbell and Perch Lundow, Entropy 21, 978 (2019).
  * (25) J. Salas, and A. D. Sokal, J. Stat. Phys. 98, 551 (2000).
  * (26) E. J. Flores-Sola, B. Berche, R. Kenna and M. Weigel, Eur. Phys. J. B. 88, 28 (2015).
  * (27) L. Jones and A. P. Young, Phys. Rev. B 71, 174438 (2005).
  * (28) N. W. Dalton and D. W. Wood, J. Math. Phys. 10, 1271 (1969).
  * (29) A. Malakis, P. Kalozoumis and N. Tyraskis, Eur. Phys. J. B 50, 63 (2006).

Figure 1: (a) Binder Cumulant; (b) and the second-order correlation length
using four lattice sizes, for fixed ferromagnetic $J_{1}$ and $J_{2}$
couplings (case $p=1$ in Eq.(2)). The scale invariance at the critical
temperature is clearly seen at the crossing point in (a) and (b). The critical
temperature is $T_{c}\approx 5.26$.

Figure 2: (a) Binder Cumulant; (b) and the second-order correlation length
using four lattice sizes, for fixed ferromagnetic $J_{1}$ and
antiferromagnetic $J_{2}$ couplings (case $p=0$ in Eq.(2)). Note that the
critical temperature is close to $T=2.08$ (see reference malakis_2006 ).

Figure 3: (a) Binder Cumulant; (b) and the second-order correlation length
using four lattice sizes, for fixed ferromagnetic $J_{1}$ couplings, and
random $J_{2}$ couplings obeying the PDF given in Eq.(2), for $p=0.6$. This is
a nonequilibrium critical steady state. The scale invariance at the critical
temperature is clearly seen at the crossing point in (a) and (b).

Figure 4: Correlation length scanned at a fixed very low temperature, for
fixed ferromagnetic $J_{1}$ couplings, and random $J_{2}$ couplings obeying
the PDF given in Eq.(2). In (a) the critical probability of the SAF-P
transition is about $p=0.078$. In (b) the critical probability of the F-P
transition is about $p=0.57$.

Figure 5: Phase diagram of the model given by the Hamiltonian given in Eq.(1),
for fixed ferromagnetic $J_{1}$ couplings, and random $J_{2}$ bonds obeying
the PDF given in Eq.(2). Critical points for $p=0$ and $p=1$ are at
equilibrium, whereas for $0<p<1$, the frontier points are at nonequilibrium
steady states.

Figure 6: Phase diagram of the model given by the Hamiltonian given in Eq.(1),
for fixed ferromagnetic $J_{1}$ couplings, and random $J_{2}$ bonds obeying
the PDF given in Eq.(5), when $\sigma=0.1$. The vertical axis is the
temperature and the horizontal one is the probability $p$, which is the weight
of the first gaussian in Eq.(3). There is only one frontier separating the F
and P phases, since any amount of $\sigma$ causes de destruction of the SAF
phase.

Figure 7:  (a) The critical probability of the phase diagram shown in Fig6
versus the parameter $\sigma$ given in Eq.(3). Note the discontinuous fall in
$p_{c}$ when $\sigma$ assumes a finite value. In (b) is shown the PDF given in
Eq.(3) for three different values of the pair ($p$, $\sigma$). Note that there
is a weak overlap of the two peaks, for $p_{c}=0.33$ and $\sigma=0.3$ causing
the first inflection point in the curve shown in (a). For $p_{c}=0.987$ and
$\sigma=2.5$ the left gaussian is just irrelevant, and the right one is very
extended.

Figure 8:  (a) The critical temperature of the phase diagram shown in Fig6
versus the parameter $\sigma$ given in Eq.(3), for $p=1$. Note the
discontinuous fall in $T_{c}$ when $\sigma$ assumes a finite value. In (b) is
shown the PDF given in Eq.(3) for three different values of $\sigma$, for
$p=1$. Note that the critical temperature reaches its maximum for $\sigma=1$,
and decreases when the left tail of the gaussian include more negative values
of $J_{2}$.

Figure 9: Phase diagram of the model given by the Hamiltonian given in Eq.(1),
for fixed antiferromagnetic $J_{1}$ interactions, and random $J_{2}$ bonds
obeying the PDF given in Eq.(5). The maxima at both frontiers are not present
in the equivalent phase diagram of Fig.5a of reference yining_2018 , which was
obtained with quenched randomness at equilibrium.

Figure 10: (a) Curves of $\xi_{SAF}/L$, where the crossing point is at
$T=2.196\pm 0.003$ (above $T_{c}(p=0)$); (b) Curves of $\xi_{AF}/L$, where the
croosing point is at $T=2.364\pm 0.004$ (above $T_{c}(p=1)$). These two
crossing points correspond to the maximum of each frontier in Fig.9.
Furthermore, in (a) $\xi_{SAF}/L$ is $1.467\pm 0.008$ at the point of
intersection , whereas in (b) $\xi_{AF}/L$ is $0.904\pm 0.003$ (which is the
Ising universal amplitude $x^{*}$ salas2000 ). Dashed lines are guides to the
eye.

Figure 11: Phase diagram for fixed antiferromagnetic $J_{1}$ bonds, and random
$J_{2}$ interactions obeying the PDF given in Eq.(6), when $\sigma=0.1$. The
horizontal axis is the probability $p$, which is the weight of the first
gaussian in Eq.(6). There is only one frontier separating the AF and P phases.
The SAF phase can not be formed for $\sigma\neq 0$.

Figure 12:  (a) The critical probability of the phase diagram shown in Fig11
versus the parameter $\sigma$ given in Eq.(6). Note the discontinuous fall in
$p_{c}$ when $\sigma$ assumes a finite value. In (b) is shown the PDF given in
Eq.(6) for three different values of the pair ($p$, $\sigma$).

Figure 13:  (a) The critical temperature of the phase diagram shown in Fig11
versus the parameter $\sigma$ given in Eq.(6), for $p=1$. Note that $\sigma$
does not affect the value of the critical temperature for $\sigma\lesssim
0.35$. In (b) is shown the PDF given in Eq.(6) for two different values of
$\sigma$, for $p=1$. Note that the tails of the gaussian for $\sigma=0.35$
include $|J_{2}|>1$, whereas for $\sigma=0.2$ they die for $|J_{2}|<1$.

Figure 14: Curves of the correlation length for different sizes of the
lattice. Curves in (a) and (b) intercept at the critical temperature of the
SAF-P transition in Figures 5 and 9, for $p=0$, respectively. In (a) the
crossing point is observed at $(2.0832\pm 0.0006,1.210\pm 0.005)$ ; In (b) it
is at $(2.082\pm 0.001,1.483\pm 0.008)$. Dashed lines are guides to the eye,
through which we may note that $x^{*}_{SAF}$ is different in (a) and (b).

Figure 15:  Extrapolation of $\xi_{SAF}/L$ to determine the amplitude
$x^{*}_{SAF}$ for $L\to\infty$ (see the intercepts at the vertical axis). It
is based on Eq.(16). The lines correspond to the limits of the interval of the
uncertainty of $T_{c}$ (read the caption of Fig.14b). Therefore, the points
between the intercepts are $1.440<x^{*}_{SAF}<1.497$.

Figure 16: In (a) the uncertainty of the crossing point gives the critical
probability $p$ of the PDF in Eq.(5), in the interval $0.30882<p_{c}<0.30898$.
In (b) the limits of the uncertainty of $x^{*}_{SAF}$ are determined by
extrapolating $\xi_{SAF}/L$ to $L\to\infty$, for each limit of $p_{c}$. So we
have $1.53<x^{*}_{SAF}<1.55$. This interval does not intersect that of Fig.15.


# Language Identification of Bengali-English Code-Mixed data using Character &
Phonetic based LSTM Models

Soumil Mandal1, Sourya Dipta Das2, Dipankar Das3

1Department of CSE, SRM University, Chennai, India
2Department of ETCE, Jadavpur University, Kolkata, India
3Department of CSE, Jadavpur University, Kolkata, India
{soumil.mandal1, dipta.math2<EMAIL_ADDRESS>
###### Abstract

Language identification of social media text still remains a challenging task
due to properties like code-mixing and inconsistent phonetic transliterations.
In this paper, we present a supervised learning approach for language
identification at the word level of low resource Bengali-English code-mixed
data taken from social media. We employ two methods of word encoding, namely
character based and root phone based to train our deep LSTM models. Utilizing
these two models we created two ensemble models using stacking and threshold
technique which gave 91.78% and 92.35% accuracies respectively on our testing
data.

##  1 Introduction

Exploration of social media and exponential growth of data redirect us to
bring various insights on language aspects. For example, a phenomenon known as
code-mixing where people from multilingual backgrounds often mix two languages
while communicating on social media. As our traditional NLP tools often fail
or give poor results on such data Bergsma et al. (2012), it becomes an
important prerequisite to identify language tags with high accuracy. Language
identification on such data at both word and sentence level still remains as a
challenging task due to a few facts such as inconsistent phonetic
transliteration, borrowing of words, spelling errors, intra sentential change
of matrix, embedded languages as well as the use of numeric and special
characters in words. Here, we discuss our system for word level language
identification of low resource Bengali-English code-mixed data taken from
Twitter, where both are typed in Roman script. Bengali is an Indo-Aryan
language of India where 8.10%
111https://en.wikipedia.org/wiki/Languages_of_India of the total population
are the $1^{st}$ language speakers and is also the official language of
Bangladesh. The original script in which Bengali is written by locals is the
Eastern Nagari Script 222https://www.omniglot.com/writing/bengali.htm. Most of
the code-mixed data was collected from Twitter 333https://twitter.com/.
Employing LSTM, two models were developed, one using character based encoding
and the other using Bengali root phones. Combining these two models, two
ensemble models were prepared using stacking and threshold. Finally, the
performance metrics like accuracy and macro-averaged F1-Score along with
confusion matrix were used to compare the results. With respect to the related
works in Sec 2, our system managed to score comparable results using a
fraction of data.

##  2 Related Work

The LID problem in code-mixed texts has already been studied by many
researchers in the aspect of social media communication utilizing various
methods involving both supervised and unsupervised techniques. Nguyen and
Doğruöz (2013) experimented with linear chain CRFs but their contextual
features are limited to bi-grams of words. Kim et al. (2014) studied the
linguistic behavior of bilingual Twitter users from Qatar, Switzerland and
Québec, and also acknowledged that code-mixing could not be studied due to the
absence of appropriate LID tools. Das and Gambäck (2014) introduced code-
mixing index to evaluate the level of mixing of the languages and made a
system for Bengali-English LID using SVM trained on n-gram with weights,
dictionary, minimum edit distance, 7 word window for context and got an
F1-Score of 76.37%. Barman et al. (2014) uses char-n-grams, capitalization,
dictionary along with a decision tree trained on word length to train a SVM
model for code-mixed LID. They also use CRF for context capturing. Jhamtani et
al. (2014) combined to classifiers into an ensemble model for Hindi-English
code-mixed LID. The first classifier used modified edit distance, word
frequency and character n-grams as features. The second classifier used output
of the first classifier for the current word, along with language tag and POS
tag of neighboring to give the final tag. Piergallini et al. (2016) used
n-grams along with capitalization as features to train a linear model for
Swahili-English word level LID and got an accuracy of 96.5%. Jaech et al.
(2016) employed a hierarchical neural model, where a CNN was trained on
char2vec and bi-LSTM on sentence level word embeddings to create a system
which got an F1-Score of 93.1% on Spanish-English data. Rijhwani et al. (2017)
demonstrates a generalized unsupervised model based on HMM for code-mixed LID
for an arbitrarily large set of languages. Singer et al. (2003) used phones
and acoustic properties for language identification in speech data. To the
best of our knowledge, none of the works have used character or phonetic
encoding in an LSTM architecture for textual language identification task.

##  3 Data Sets

We collected transliterated Bengali words in Roman and English words from ICON
16 444http://ltrc.iiit.ac.in/icon2016/, ICON 17
555https://ltrc.iiit.ac.in/icon2017/ contests and code-mixed data used in
Mandal and Das (2018). Some additional English words were collected from
online 666https://github.com/first20hours/google-10000-english resources. The
datasets had no intersection. All words were converted to lowercase, and words
having length less than 3 or any special or numeric characters (e.g. ri8,
2toh) were removed. Also, words with $>$ 2 consecutive identical characters
was normalized to 2 consecutive identical characters. Words which had both EN
and BN components, for e.g. journey(ta), which is an EN word with BN suffix
were discarded as well. The distribution for training, development and testing
was 6632x2, 300x2 and 700x2, respectively. It should be noted that the amount
of data is comparatively less as Bengali falls in the low resource language
category when dealing with data from social media.

##  4 Support Vector Machine Model

In order to start with (baseline model), we decided to make a supervised model
using Support Vector Machine (SVM) with bag-of-words (BOW) principle. The
training and development set were merged in order to create the training data
in this case. Then character level bigrams, trigrams and quadgrams were
extracted and SVM with linear kernel implemented in the scikit-learn
777http://scikit-learn.org/stable/ package was employed for training purpose.
This system achieved an accuracy of 83.64% on the test data.

##  5 Neural Model

We moved towards implementing neural models to achieve better accuracy. Two
encoding methods for neural model were tried, namely character based and phone
based. The encoding methods we used are described in detail below.

###  5.1 Character Encoding

We decided to make a neural model based on character encoding for acquiring
patterns at an elementary level. Similar encoding methods have been tested for
character level NMT before, used in a seq2seq architecture Lee et al. (2016).
In order to build the encoder, a dictionary based on character-index was made
where index is the $n^{th}$ English alphabet. The algorithm fetches each
character at a time from the word and replaces it with the respective index,
for example, the encoding of good is represented as [7,15,15,4], for bad it
will be [2,1,4].

###  5.2 Phonetic Encoding

It is evident that different languages have different phonetic properties
888https://en.wikipedia.org/wiki/Phonetics, especially when they have distant
origins. In order to exploit this, we developed a language identification
neural model based on phonetic encoding. One of our major aims was to obtain
such encodings in terms of Bengali root phones. Thus, a phonetic library that
contains two parts, namely root phones and similar phone groups was prepared
by the authors in 3 steps. First, character level bigrams and trigrams were
collected from the Bengali words present in the training set and were sorted
in non-increasing order of frequency. We did not extend to quadgrams as we
could see by observing the data set that the users generally tend to use
maximum of 3 characters to phonetically represent the corresponding Bengali
character. Then, from the code chart of Bengali unicode characters
999https://unicode.org/charts/PDF/Unicode-10.0/, the pronunciation of each of
the Bengali characters written in Roman script was gathered and finally, a
root phone list (RP) was prepared. By referring to the (2,3)-gram frequency
table, we noticed that some of the transliterations collected from the unicode
chart are rarely used in real life for phonetic typing. Thus, such instances
were discarded from the list (e.g. nga, nya, ddha). Next, with respect to the
root phones, all the corresponding (phonetically similar) bigrams and trigrams
with high frequency were gathered and formed into groups, i.e similar phone
groups (SPG). Though some of the bigrams and trigrams had quite high
frequency, yet they were not included in SPG as they didn’t correspond to any
of the root phones, example (e.g. bhi). The phonetic library was finally
checked and verified by a linguist. Complete model consisting of RP and SPG is
shown in Table 1. In each of the phonetic groups, we had set the first phone
as the respective root phone for ease of searching and retrieval of root
index.

Root Phones
---
| aa, i, u, r, e, ai, o, au, ka, kha, ga, gha, ca, cha,
---
ja, jha, ta, tha, da, dha, na, pa, pha, ba, bha, ma,
ya, ra, la, sa, ha
Phonetic Groups
(aa, a), (i, ee), (u, w), (r, ri), (e), (ai, oi), (o, oo)
(au, ou, ow), (ka, k), (kha, kh), (ga, g), (gha, gh)
(ca, c), (cha, ch), (sa, s, sh), (jha, jh), (bha, bh, v)
(ta, t), (tha, th), (da, d), (dha, dh), (na, n), (pa, p)
(pha, ph, f), (ba, b), (ma, m), (ya, y), (ra, rh)
(la, l), (ja, j, z), (ha, h)
Table 1: Phonetic library.

The routine for encoding is based on variable length sliding window technique,
where n$\leftarrow$[3,1]. The algorithm in each of iteration takes n
characters, with n starting as 3. If the character cluster appear in phonetic
groups, we append the root phone index, else we reduces n by 1 and repeat the
same process till n = 1. If the n-gram is not found in phonetic groups, 35 (as
number of root phones is 31) is appended. From the next iteration, the window
starts from n + 1, where n denotes the ending point of the window in the last
iteration. The algorithm stops when window traverses the whole word. Pseudo-
code given below.

while index $<$ len(word)
if word[index:index+j](j$\rightarrow$3,2,1) in SPG
enc.append(root_phone_index), index += j
else enc.append(35), index += 1
return enc

On encoding the words in the combined dataset using the above algorithm, the
number of times 35 occurred (i.e n-gram absent in the phonelib is seen) was
found to be 82 and 1 for EN and BN data respectively. Some examples of word
encodings are, khabar meaning ’food’ will be [10,24,4], which is essentially
kha-ba-r. Also khbr, which is the same word but is less seen in terms of
phonetic typing, produces identical encoding. Though this transforming ability
performs quite well on many of such instances, it fails on some instances as
well like the word korchi produces [9,7,4,14,2] while krci results in
[9,4,13,2] even though they are the phonetic transliterations of the same
word, where the latter is less accurately typed. The normalized root phones
frequency graph is shown in Fig 1. There, y-axis denotes the normalized
frequency while the x-axis denotes the root phones (same order as in Table 1).
From the table we can see 3 spikes in the starting, which are essentially the
frequency of the phones i,e,o. Thus we can draw the inference, that these 3 RP
are quite essential while transliterating BN in Roman. From the rest of the
graph, we can see a clear distinction between the blue and green line,
implying how EN and BN is different in terms of the root phones.

Figure 1: Normalized root phones frequency.

##  6 Training

We adopted the Long Short Term Memory version of the RNN architecture
101010https://keras.io/ to train both of our models. We specifically chose
this because it outperforms almost all the other models in sequence learning
Greff et al. (2017). Two LSTM models were built, one for each type of
encoding. We set the same parameters for both, with sigmoid activation
function and adam optimizer whereas the loss function used was binary cross-
entropy, number of epochs was set to 500 and batch size was set to 1658. All
the other parameters were kept at default. The significant difference between
the two models was the architecture. For character encoding, the layer sizes
were 15-35-25-1 and for phonetic it was 15-15-40-1. The numerics indicate here
that the input and the output layer has 15 and 1 node(s) respectively, and the
other two are the number of nodes in two hidden layers. The encodings were
padded to length 15 in order to make them uniform and was finally converted
into one-hot vectors. The target label for Bengali and English was set to 0
and 1, respectively, thus the output of the trained models was fuzzy values
between them (this can be visualized clearly in the scatter plots given in the
next section).

##  7 Evaluation

For both the character and phonetic models, first a simple roundup logic was
used for testing on the development data (dev_round). As from the fuzzy value
scatter plots we could see that theres a scope of improvement, we decided to
use brute force for finding a threshold on the development data (dev_thresh)
in order to tune the trained model. Finally using this threshold, we evaluated
the model on our testing data (test_thresh). The details of both the models
along with the ensemble model is given below.

###  7.1 Character Model

The results of the character model are shown in Table 2 and the fuzzy outputs
scatter plot on the development data is shown in Fig 2. Simple round up on the
fuzzy output gave an accuracy of 91.50%. Using brute force, the threshold
where the accuracy is peaking on the development data was found to be at
$\theta$ $\leq$ 0.92. This showed a slight improvement in accuracy by 0.66%.
Finally, on test data, this tuned model gave an accuracy of 91.71%.

Figure 2: Scatter plot of character model fuzzy values on development data. Model | Acc | Prec | Rec | F1
---|---|---|---|---
dev_round | 91.50 | 87.84 | 96.33 | 91.48
dev_thresh | 92.16 | 91.47 | 93.00 | 92.16
test_thresh | 91.71 | 91.59 | 91.85 | 91.71
Table 2: Character model results (in %).

###  7.2 Phonetic Model

Results of the phonetic model are shown in Table 3. On simple roundup on
development data, the accuracy was not very good (82.5%), but using a similar
threshold technique as used in the character model, the accuracy improved
significantly (6.16%) and reached 88.66%. The threshold in this case was
calculated to be $\theta$ $\leq$ 0.95. The scatter plot of the fuzzy values is
shown in Fig 3. Using threshold, the accuracy achieved on the test data was
90.42% and a precision of 91.74% which is slightly better than the character
model.

Figure 3: Scatter plot of phonetic model fuzzy values on development data. Model | Acc | Prec | Rec | F1
---|---|---|---|---
dev_round | 82.50 | 75.06 | 97.33 | 82.10
dev_thresh | 88.66 | 87.66 | 90.00 | 88.66
test_thresh | 90.42 | 91.74 | 88.85 | 90.42
Table 3: Phonetic model results (in %).

###  7.3 Ensemble Model

Two ensemble methods were tried, using stacking and thresholding. For the
former, logistic regression was used as the combiner algorithm and was trained
on the fuzzy values given by the two models on the development data. For
threshold, mean of the two fuzzy values was taken and by using a similar brute
force technique as used in char and phonetic model, the threshold was found to
be $\theta$ $\leq$ 0.9. The architecture is shown below in Fig 4.

Figure 4: Ensemble architecture.

The results compared with the baseline model is shown in Table 4 and the
scatter plot of means is shown in Fig 5.

Figure 5: Scatter plot of ensemble model fuzzy values on development data. Model | Acc | Prec | Rec | F1
---|---|---|---|---
svm_baseline | 83.64 | 83.88 | 83.28 | 83.64
ensem_stack | 91.78 | 89.58 | 94.57 | 91.77
ensem_thresh | 92.35 | 94.99 | 89.42 | 92.35
Table 4: Ensemble model results (in %).

The stacking method showed slight improvement and the accuracy achieved was
91.78%. On the other hand, accuracy improved by 0.57% using threshold
technique. Also, we can see the improvement using the neural networks is quite
noticeable, i.e an increase in accuracy by 8.71% from the SVM n-gram based
model.

###  7.4 Error Analysis

We started our error analysis by preparing confusion matrices (CM) of the four
models. The values observed on test data are shown in Table 5. Here, BN
correctly predicted was considered as TP and EN was considered as TN.

Confusion Matrices
---
1 | BN |  EN | 2 | BN | EN
BN | 641 | 59  | BN | 644 | 56
EN | 57 | 643  | EN | 78 | 622
3 | BN |  EN | 4 | BN | EN
BN | 623 | 77  | BN | 667 | 33
EN | 38 | 662  | EN | 74 | 626
Table 5: Confusion matrices. (1-character, 2-phonetic, 3-ensem_stack,
4-ensem_thresh)

From CM 1 given in Table 5, we found that the character model is quite
balanced in terms of biasness as the values of TP and TN are quite close.
Unlike the character model, we can see from CM 2 (phonetic model) that TP is
much more than TN, thus the precision is higher by 0.15% compared to the
character model. From CM 3, we observed that TN is much higher compared to TP
(by 39) and thus, it can be inferred that the logistic regression model is
biased towards EN. On the contrary, from CM 4, we can claim that TP is much
higher compared to TN (by 41). This is due to the fact that the regression
range of the fuzzy outputs from the models for BN is much higher compared to
that of EN, and thus calculating threshold keeping accuracy as the metric
favors TP.

##  8 Observations

One of the primary drawbacks of our system is that it is based on word level
and not sentence, thus, it fails in capturing context information. This
problem is quite evident when words with similar spellings but belonging to
different languages are assigned with same tag. E.g, in our test data, the
word ”choke” was present in both EN and BN (meaning eyes) data but was always
tagged as BN. The system is not designed to handle words with numeric or
special characters (e.g. ri8 instead of right), and also elongated words (e.g.
goood instead of good) though a simple normalization prior to feeding would
solve the issue.

##  9 Conclusion & Future Work

In the present work, we have build a system for word level LID. A new method
for word encoding has been introduced using root phones. Employing character
and phonetic encoding methods, we designed two deep LSTM models. Finally, two
ensemble models were designed using stacking and threshold techniques. The
stacking method achieved an accuracy of 91.78% while on the other hand the
threshold model obtained the best accuracy with a score of 92.35%. Considering
that our training data was quite low $\approx$13k, the architecture seems to
be quite promising. In future, our aim is to gather more data from social
media. This would not only help in the generic improvement of accuracy, but
the n-gram statistics may also be useful to enrich and fine tune our phonetic
library even more, thus improving the encoding quality. Also, we would like to
test a similar methodology on other language pairs as well to evaluate its
applicability.

## References

  * Barman et al. (2014) Utsab Barman, Amitava Das, Joachim Wagner, and Jennifer Foster. 2014.  Code mixing: A challenge for language identification in the language of social media.  In Proceedings of the first workshop on computational approaches to code switching. pages 13–23.
  * Bergsma et al. (2012) Shane Bergsma, Paul McNamee, Mossaab Bagdouri, Clayton Fink, and Theresa Wilson. 2012.  Language identification for creating language-specific twitter collections.  In Proceedings of the second workshop on language in social media. Association for Computational Linguistics, pages 65–74.
  * Das and Gambäck (2014) Amitava Das and Björn Gambäck. 2014.  Identifying languages at the word level in code-mixed indian social media text .
  * Greff et al. (2017) Klaus Greff, Rupesh K Srivastava, Jan Koutník, Bas R Steunebrink, and Jürgen Schmidhuber. 2017.  Lstm: A search space odyssey.  IEEE transactions on neural networks and learning systems 28(10):2222–2232.
  * Jaech et al. (2016) Aaron Jaech, George Mulcaire, Mari Ostendorf, and Noah A Smith. 2016.  A neural model for language identification in code-switched tweets.  In Proceedings of The Second Workshop on Computational Approaches to Code Switching. pages 60–64.
  * Jhamtani et al. (2014) Harsh Jhamtani, Suleep Kumar Bhogi, and Vaskar Raychoudhury. 2014.  Word-level language identification in bi-lingual code-switched texts.  In Proceedings of the 28th Pacific Asia Conference on Language, Information and Computing.
  * Kim et al. (2014) Suin Kim, Ingmar Weber, Li Wei, and Alice Oh. 2014.  Sociolinguistic analysis of twitter in multilingual societies.  In Proceedings of the 25th ACM conference on Hypertext and social media. ACM, pages 243–248.
  * Lee et al. (2016) Jason Lee, Kyunghyun Cho, and Thomas Hofmann. 2016.  Fully character-level neural machine translation without explicit segmentation.  arXiv preprint arXiv:1610.03017 .
  * Mandal and Das (2018) Soumil Mandal and Dipankar Das. 2018.  Analyzing roles of classifiers and code-mixed factors for sentiment identification.  arXiv preprint arXiv:1801.02581 .
  * Nguyen and Doğruöz (2013) Dong Nguyen and A Seza Doğruöz. 2013.  Word level language identification in online multilingual communication.  In Proceedings of the 2013 Conference on Empirical Methods in Natural Language Processing. pages 857–862.
  * Piergallini et al. (2016) Mario Piergallini, Rouzbeh Shirvani, Gauri S Gautam, and Mohamed Chouikha. 2016\.  Word-level language identification and predicting codeswitching points in swahili-english language data.  In Proceedings of the Second Workshop on Computational Approaches to Code Switching. pages 21–29.
  * Rijhwani et al. (2017) Shruti Rijhwani, Royal Sequiera, Monojit Choudhury, Kalika Bali, and Chandra Shekhar Maddila. 2017.  Estimating code-switching on twitter with a novel generalized word-level language detection technique.  In Proceedings of the 55th Annual Meeting of the Association for Computational Linguistics (Volume 1: Long Papers). volume 1, pages 1971–1982.
  * Singer et al. (2003) Elliot Singer, Pedro A Torres-Carrasquillo, Terry P Gleason, William M Campbell, and Douglas A Reynolds. 2003.  Acoustic, phonetic, and discriminative approaches to automatic language identification.  In Eighth European conference on speech communication and technology.


# Difference independence of the Euler gamma function

Qiongyan Wang1, Xiao Yao2 1 School of Mathematical sciences, Peking
University, Beijing, 100871, P.R. China 2 School of Mathematical sciences and
LPMC, Nankai University, Tianjin, 300071, P.R. China
<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract.

In this paper, we established a sharp version of the difference analogue of
the celebrated Hölder’s theorem concerning the differential independence of
the Euler gamma function $\Gamma$. More precisely, if $P$ is a polynomial of
$n+1$ variables in $\mathbb{C}[X,Y_{0},\dots,Y_{n-1}]$ such that

$P(s,\Gamma(s+a_{0}),\dots,\Gamma(s+a_{n-1}))\equiv 0$

for some $(a_{0},\dots,a_{n-1})\in\mathbb{C}^{n}$ and
$a_{i}-a_{j}\notin\mathbb{Z}$ for any $0\leq i<j\leq n-1$, then we have

$P\equiv 0.$

Our result complements a classical result of algebraic differential
independence of the Euler gamma function proved by Hölder in 1886, and also a
result of algebraic difference independence of the Riemann zeta function
proved by Chiang and Feng in 2006.

Corresponding author: Xiao Yao

2010 Mathematics Subject Classification. 11M06, 39A05.

The research was partially supported by National Key R&D Program of China
(2020YFA0713300) and NSFC of China(No.11901311).

Keywords. Algebraic difference independence; Euler gamma function; Algebraic
difference equations.

##  1\. Introduction

A classical theorem of Hölder [8] states that the Euler gamma function

$\Gamma(s)=\int_{0}^{+\infty}t^{s-1}e^{-t}dt,~{}~{}~{}\Re s>0$

which can be analytically continued to the whole complex plane $\mathbb{C}$,
does not satisfy any nontrivial algebraic differential equation whose
coefficients are polynomials in $\mathbb{C}$. We state it in the following.

######  Theorem A.

Let $P$ be a polynomial of $n+1$ variables in
$\mathbb{C}[X,Y_{0},\dots,Y_{n-1}]$. Assume that

$P(s,\Gamma(s),\dots,\Gamma^{(n-1)}(s))\equiv 0,$

then we have

$P\equiv 0.$

To the best of our knowledge, the Euler gamma function $\Gamma$ seems to be
the first known example which satisfies the algebraic differential
independence property in the literature. It is well known that the Riemann
zeta function $\zeta$ is associated with $\Gamma$ by the famous Riemann
functional equation

(1) $\displaystyle\zeta(1-s)=2^{1-s}\pi^{-s}\cos\frac{\pi
s}{2}\Gamma(s)\zeta(s).$

Motivated by the Riemann functional equation, it is natural to consider the
algebraic differential independence property for the Riemann zeta function.
The study of the algebraic differential independence of the Riemann zeta
function $\zeta$ can be dated back to Hilbert. In [7], he conjectured that
Hölder’s result can be extended to the Riemann zeta function $\zeta$. Later,
this conjecture was verified by Ostrowski in [15].

Bank and Kaufman [2, 3] made the following celebrated generalizations of
Hölder’s result.

######  Theorem B.

Let $P$ be a polynomial in $K[X,Y_{0},\dots,Y_{n-1}]$, where $K$ is the field
of all meromorphic functions such that the Nevanlinna’s characteristic
$T(r,f)=o(r)$ as $r$ goes to infinity for any $f$ in $K$. Assume that

$P(s,\Gamma(s),\dots,\Gamma^{(n-1)}(s))\equiv 0,$

then we have

$P\equiv 0.$

For the Nevanlinna characteristic $T(r,f)$, we refer to Hayman’s book [6] for
a detailed introduction. Since $\Gamma$ and $\zeta$ appeared very naturally in
Riemann functional equation (1), Markus in [14] posted an open problem to
study the joint algebraic differential independence of $\Gamma$ and $\zeta$.
We refer the readers to the references [10, 11, 12, 13] for the recent
developments in this direction.

It is interesting to study the algebraic difference independence of $\zeta$ or
$\Gamma$. Feng and Chiang proved the following result.

######  Theorem C.

Let $P$ be a polynomial of $n+1$ variables in
$\mathbb{C}[X,Y_{0},\dots,Y_{n-1}]$ and $s_{0},\dots,s_{n-1}$ be $n$ distinct
numbers in $\mathbb{C}$. Assume that

$P(s,\zeta(s+s_{0}),\dots,\zeta(s+s_{n-1}))\equiv 0,$

then we have

$P\equiv 0.$

Chiang and Feng’s result extended a result of Ostrowski in [15] where the
assumption of $s_{0},\dots s_{n-1}$ are $n$ distinct real numbers are needed.
Indeed, Chiang and Feng proved that Theorem C also holds under the same
assumption in Theorem B , we refer the interested readers to [4] for the
details. Here, we also mention two remarkable universality results due to
Voronin in 1970s for the differential case [17] and the difference case [18].
We refer to [16] for the detailed introduction of the recent developments in
this direction.

To the best of our knowledge, the topic of the algebraic difference
independence of the Euler gamma function was first addressed by Hardouin in
[5] in the framework of difference Galois theory. Motivated by the
multiplication theorem of Euler gamma function

(2)
$\Gamma(ns)=n^{ns-\frac{1}{2}}(2\pi)^{\frac{1-n}{2}}\prod_{j=0}^{n-1}\Gamma(s+\frac{j}{n}),$

Hardouin proved the following result.

######  Theorem D ([5]).

Let $a_{0},\dots,a_{n-1}$ be $n$ complex numbers in $\mathbb{C}$,
$b_{0},\dots,b_{m-1}(\geq 2)$ be $m$ integers such that $\\{a_{j}(mod\
1)\\}_{j=0}^{n-1}$ and $\\{\sum\limits_{l=0}^{b_{j}-1}\frac{l}{b_{j}}(mod\
1)\\}_{j=0}^{m-1}$ are $\mathbb{Z}$-linearly independent. Assume that

$P(s,\Gamma(s+a_{0}),\dots,\Gamma(s+a_{n-1}),\Gamma(b_{0}s),\dots,\Gamma(b_{m-1}s))\equiv
0$

for some polynomial $P$, then we have

$P\equiv 0.$

Hardouin’s proof relies on Kolchin’s type theorem in an essential way. See
also in [1] for a detailed disccusion of Kolchin’s type theorem and several
powerful applications in algebraic independence problems.

Our starting point is another well known difference equation of $\Gamma$,

(3) $\Gamma(s+1)=s\Gamma(s).$

This may be the obvious obstruction for us to study the algebraic difference
independence of the Euler gamma function $\Gamma$. One can not expect to
obtain Theorem B for $\Gamma$ directly. While in this paper, we will show that
the machinery exhibited in (3) is the only obstruction to get the algebraic
difference independence of $\Gamma$. Now, we state our main result in the
following. In this paper, we will use an elementary method inspired by [8, 15]
to prove our main result, which avoids the advanced difference Galois theory.
This may be of independent interest.

We define

(4)
$\mathcal{H}:=\\{(a_{0},\dots,a_{n-1})\in\mathbb{C}^{n}:a_{i}-a_{j}\notin\mathbb{Z}\
\text{for any}\ 0\leq i<j\leq n-1\\}.$

Now, we state our main result in the following.

######  Theorem 1.

Let $P$ be a polynomial of $n+1$ variables in
$\mathbb{C}[X,Y_{0},\dots,Y_{n-1}]$. Assume that

$P(s,\Gamma(s+a_{0}),\dots,\Gamma(s+a_{n-1}))\equiv 0$

for some $(a_{0},\dots,a_{n-1})\in\mathcal{H}$, then we have

$P\equiv 0.$

We remark that we can also use Theorem D to recover part of the result of
Theorem 1 under the same condition of $(a_{j})_{j=0}^{n-1}$ and also $m=0$ in
Theorem D. While, it can not completely recover Theorem 1, since the condition
in Theorem 1 is sharp. Our result complements the classical result of
algebraic differential independence of Euler gamma function proved by Hölder
[8] in 1886, and also a result of algebraic difference independence of Riemann
zeta function proved by Chiang and Feng [4] in 2006.

######  Corollary 1.

Let $P$ be a polynomial of $n+1$ variables in
$\mathbb{C}[X,Y_{0},\dots,Y_{n-1}]$. Assume that

$P(s,\Gamma(s),\dots,\Gamma(s+(n-1)\alpha))\equiv 0$

for some $\alpha\not\in\mathbb{Q}$, then we have

$P\equiv 0.$

######  Remark 1.

Theorem 1 can be seen as a difference version of the Hölder’s theorem. The
identity (3) shows that the discussion restricted to $\mathcal{H}$ is
necessary.

We can also extend Theorem 1 to the setting of $K[X,Y_{0},\dots,Y_{n-1}]$
where $K$ is the field of all meromorphic functions such that the Nevanlinna’s
characteristic $T(r,f)=o(r)$ as $r$ goes to infinity for any $f$ in $K$.
While, we will not address it in this paper.

By Theorem 1 and the Euclidean’s algorithm, it is not hard to give the
following two examples.

######  Example 1.

Let $P=P(X,Y,Z)$ be a polynomial of $3$ variables in $\mathbb{C}[X,Y,Z]$.
Assume that

$P(s,\Gamma(s+a_{0}),\Gamma(s+a_{1}))\equiv 0,$

then

$P\equiv 0,$

unless $a_{1}-a_{0}\in\mathbb{Z}$. In the latter case, if $\Re a_{0}<\Re
a_{1}$, $P$ can be divided by the polynomial
$R(X,Y,Z)=Z-(X+a_{0})\dots(X+a_{1}-1)Y$.

######  Example 2.

Let $P(X,Y,Z,W)=YW-Z^{2}-YZ$ in $\mathbb{C}[X,Y,Z,W]$. We have

$P(s,\Gamma(s),\Gamma(s+1),\Gamma(s+2))\equiv 0.$

$P$ belongs to the ideal

$<W-(X+1)Z,Z-XY>$

generated by $W-(X+1)Z$ and $Z-XY$ in $\mathbb{C}[X,Y,Z,W]$. Furthermore, $P$
can be written by

$P(X,Y,Z,W)=Y(W-(X+1)Z)+Z(XY-Z).$

######  Remark 2.

Indeed, inspired by Example 1 and Example 2, we can apply Theorem 1 and the
Euclidean’s algorithm again to give a complete characterization of the
following set

$\mathcal{I}:=\\{P\in\mathbb{C}[X,Y_{0},\dots,Y_{n-1}]:P(s,\Gamma(s+a_{0}),\dots,\Gamma(s+a_{n-1}))\equiv
0\\}$

without any assumption on $a_{0},\dots,a_{n-1}$. While, we will not discuss it
in this paper.

##  2\. Proof of Theorem 1

In order to prove Theorem 1, we need introduce a lexicographic order between
any two monomials $Y_{0}^{i_{0}}\dots Y_{n-1}^{i_{n-1}}$ and
$Y_{0}^{j_{0}}\dots Y_{n-1}^{j_{n-1}}$ in $\mathbb{C}[Y_{0},\dots,Y_{n-1}]$,
which plays an important role in our proof. And this strategty was inspired by
Ostrowski’s proof of Hölder’s classical proof in [15]. It also shares some
spirit of Kolchin’s type theorem which was used in [5, 1]

We first introduce an order for the $n$ symbols $Y_{0},\dots,Y_{n-1}$,

(5) $Y_{0}\prec Y_{1}\prec\dots\prec Y_{n-1}.$

This can be used to induce a lexicographic order between any two monomials
$Y_{0}^{i_{0}}\dots Y_{n-1}^{i_{n-1}}$ and $Y_{0}^{j_{0}}\dots
Y_{n-1}^{j_{n-1}}$. We still denote it by $\prec$ to simplify the notation. We
define it in the following,

  1. case 1:

$Y_{0}^{i_{0}}\dots Y_{n-1}^{i_{n-1}}=Y_{0}^{j_{0}}\dots Y_{n-1}^{j_{n-1}}$ if
$i_{k}=j_{k}$ for $k=0,\dots,n-1$;

  2. case 2:

$Y_{0}^{i_{0}}\dots Y_{n-1}^{i_{n-1}}\prec Y_{0}^{j_{0}}\dots
Y_{n-1}^{j_{n-1}}$ if $i_{0}<j_{0}$ or there exists $1\leq k\leq n-1$ such
that

$i_{0}=j_{0},\dots,i_{k-1}=j_{k-1},i_{k}<j_{k};$

  3. case 3:

$Y_{0}^{j_{0}}\dots Y_{n-1}^{j_{n-1}}\prec Y_{0}^{i_{0}}\dots
Y_{n-1}^{i_{n-1}}$ can be defined similarly as in case 2.

For any nonzero polynomial $P=P(X,Y_{0},\dots,Y_{n-1})$ in
$\mathbb{C}[X,Y_{0},\dots,Y_{n-1}]$, we write it by

(6) $P=\sum_{i=(i_{0},\dots,i_{n-1})}\Phi_{i}(X)Y_{0}^{i_{0}}\dots
Y_{n-1}^{i_{n-1}},$

where $\Phi_{i}(X)\in\mathbb{C}[X]$ and $\Phi_{i}(X)\neq 0$. The highest term
of $P$ is defined by the maximal element in $\mathcal{T}_{P}$ with respect to
the lexicographic order $\prec$ introduced above, where

(7) $\mathcal{T}_{P}:=\\{Y_{0}^{i_{0}}\dots Y_{n-1}^{i_{n-1}}:\Phi_{i}(X)\
Y_{0}^{i_{0}}\dots Y_{n-1}^{i_{n-1}}\ \text{appeared in}\ \eqref{eq-
sum-P}\\}.$

For any monomial $L=Y_{0}^{i_{0}}Y_{1}^{i_{1}}\dots Y_{n-1}^{i_{n-1}}$, we
define its degree $\deg(L)$ by

$\deg(L):=\sum_{k=0}^{n-1}i_{k}.$

The height of $P$ is defined by the degree of the highest term of $P$.

Now, we will prove Theorem 1.

###### Proof.

Let

(8)
$\mathcal{S}:=\\{P\in\mathbb{C}[X,Y_{0},\dots,Y_{n-1}]:P(s,\Gamma(s+a_{0}),\dots,\Gamma(s+a_{n-1}))\equiv
0\\}.$

We will prove Theorem 1 by contradiction. We assume that
$\mathcal{S}\neq\\{0\\}$. By our assumption, there exists a nonzero polynomial

$Q=\sum\limits_{i=(i_{0},\dots,i_{n-1})}\Psi_{i}(X)Y_{0}^{i_{0}}\dots
Y_{n-1}^{i_{n-1}},$

which is of the lowest height in $\mathcal{S}\backslash\\{0\\}$ with
$\Psi_{j}(X)Y_{0}^{j_{0}}\dots Y_{n-1}^{j_{n-1}}$ being its highest term for
some $j=(j_{0},\dots,j_{n-1})$. Moreover, we also make the following
assumption.

Assumption LD: The nonzero polynomial $\Psi_{j}(X)$ appearing in the highest
term of $Q$ is also of the lowest degree.

Let

(9)
$T(X,Y_{0},\dots,Y_{n-1}):=Q(X+1,(X+a_{0})Y_{0},\dots,(X+a_{n-1})Y_{n-1}).$

Noting that

$Q(s,\Gamma(s+a_{0}),\dots,\Gamma(s+a_{n-1}))\equiv 0,$

we have

$T(s,\Gamma(s+a_{0}),\dots,\Gamma(s+a_{n-1}))\equiv 0$

by (3). And the highest term of $T$ is $\hat{\Psi}_{j}(X)Y_{0}^{j_{0}}\dots
Y_{n-1}^{j_{n-1}}$, where

$\hat{\Psi}_{j}(X):=\Psi_{j}(X+1)(X+a_{0})^{j_{0}}\dots(X+a_{n-1})^{j_{n-1}}.$

It follows from the Euclidean’s algorithm, there exist two polynomials
$R=R(X)$ and $U=U(X)$ in $\mathbb{C}[X]$ such that

$\hat{\Psi}_{j}=R\Psi_{j}+U,$

where either $U=0$ or $0<\deg U<\deg\Psi_{j}$. It is easy to see that $\deg
R\geq 1$.

We claim that $U=0$. Otherwise, we know that the polynomial

$H(X,Y_{0},\dots,Y_{n-1}):=T(X,Y_{0},\dots,Y_{n-1})-R(X)Q(X,Y_{0},\dots,Y_{n-1})$

is in $\mathcal{S}$. It follows that the highest term of $H$ is

$U(X)Y_{0}^{j_{0}}\dots Y_{n-1}^{j_{n-1}}$

and $0<\deg U<\deg\Psi_{j}$. Thus, $H\neq 0$, which contradicts the choice of
$Q$ and Assumption LD. Now, we have $U=0$.

Since $U=0$, we see that the highest term of $H$ is less than the highest term
of $Q$ if $H\neq 0$. This again contradicts our choice of $Q$. Thus, we get
$H=0$. That is,

(10) $T(X,Y_{0},\dots,Y_{n-1})=R(X)Q(X,Y_{0},\dots,Y_{n-1}).$

We first assume that there exists $\beta\notin\Lambda:=\\{-a_{k}:0\leq k\leq
n-1\\}$ such that $R(\beta)=0$. By (9) and (10), we get

$Q(\beta+1,(\beta+a_{0})Y_{0},\dots,(\beta+a_{n-1})Y_{n-1})=0$

in $\mathbb{C}[Y_{0},\dots,Y_{n-1}]$. This implies that

$Q(\beta+1,Y_{0},\dots,Y_{n-1})=\sum\limits_{i=(i_{0},\dots,i_{n-1})}\Psi_{i}(\beta+1)Y_{0}^{i_{0}}\dots
Y_{n-1}^{i_{n-1}}=0$

in $\mathbb{C}[Y_{0},\dots,Y_{n-1}]$. Thus, we have

$\Psi_{i}(\beta+1)=0$

for all $i$, which implies that each $\Psi_{i}(X)$ can be divided by
$X-\beta-1$. This contradicts our assumption that $\Psi_{j}$ is of the lowest
degree.

Hence, each root of $R$ lies in $\Lambda$. Without loss of generality, we
assume that $R(-a_{0})=0$. Thus, we get

$Q(-a_{0}+1,0,(a_{1}-a_{0})Y_{1},\dots,(a_{n-1}-a_{0})Y_{n-1})=0$

by (9) and (10). Recalling that $a_{j}-a_{0}\notin\mathbb{Z}$ for any $j\neq
0$, we have

(11) $Q(-a_{0}+1,0,Y_{1},\dots,Y_{n-1})=0.$

Taking $X=-a_{0}+1$, $Y_{0}=0$ in (9) and (10), we get

$\displaystyle
Q(-a_{0}+2,0,(a_{1}-a_{0}+1)Y_{1},\dots,(a_{n-1}-a_{0}+1)Y_{n-1})$
$\displaystyle=$ $\displaystyle
R(-a_{0}+1)Q(-a_{0}+1,0,Y_{1},\dots,Y_{n-1})=0$

by (11). Noting that $a_{j}-a_{0}\notin\mathbb{Z}$ for any $j\neq 0$ again, we
obtain

$Q(-a_{0}+2,0,Y_{1},\dots,Y_{n-1})=0$

in $\mathbb{C}[Y_{0},\dots,Y_{n-1}]$. By induction, we can prove that for any
$m\in\mathbb{N}$,

$Q(-a_{0}+m,0,Y_{1},\dots,Y_{n-1})=0$

in $\mathbb{C}[Y_{0},\dots,Y_{n-1}]$. It follows by the fundamental theorem of
algebra, we get

$Q(X,0,Y_{1},\dots,Y_{n-1})=0$

in $\mathbb{C}[X,Y_{0},\dots,Y_{n-1}]$. Thus, we proved that $Q$ can be
divided by the monomial $Y_{0}$, which contradicts the assumption that $Q$ is
of the lowest height in $\mathcal{S}$.

Now, we finish the proof of Theorem 1. ∎

## References

  * [1] B. Adamczewski, J. P. Bell and E. Delaygue. Algebraic independence of $G$-functions and congruences à la Lucas. Ann. Sci. Éc. Norm. Supér., 52(4), 515–559, 2019.
  * [2] S. B. Bank and R. P. Kaufman.  A note on Hölder’s theorem concerning the gamma function.  Math. Ann., 232(2):115–120, 1978.
  * [3] S. B. Bank and R. P. Kaufman.  On differential equations and functional equations.  J. Reine Angew. Math., 311(312):31–41, 1979.
  * [4] Y. M. Chiang and S. J. Feng.  Difference independence of the Riemann zeta function.  Acta Arith., 125(4):317–329, 2006.
  * [5] C. Hardouin. Hypertranscendance des systèmes aux différences diagonaux. Compos. Math. 144(3), 565–581, 2008.
  * [6] W. K. Hayman.  Meromorphic functions.  Oxford Mathematical Monographs. Clarendon Press, Oxford, 1964.
  * [7] D. Hilbert.  Mathematical problems.  Bull. Amer. Math. Soc., 8(10):437–479, 1902.
  * [8] O. Hölder.  Über die eigenschaft der $\gamma$-funktion, keiner algebraischen differentialgleichung zu genügen.  Math. Ann., 28(10):1–13, 1886.
  * [9] I. Laine.  Nevanlinna theory and complex differential equations, volume 15 of De Gruyter Studies in Mathematics.  Walter de Gruyter & Co., Berlin, 1993.
  * [10] B. Q. Li and Z. Ye.  On differential independence of the Riemann zeta function and the Euler gamma function.  Acta Arith., 135(4):333–337, 2008.
  * [11] B. Q. Li and Z. Ye.  Algebraic differential equations concerning the Riemann zeta function and the Euler gamma function.  Indiana Univ. Math. J., 59(4):1405–1415, 2010.
  * [12] B. Q. Li and Z. Ye.  Algebraic differential equations with functional coefficients concerning $\zeta$ and $\Gamma$.  J. Differential Equations, 260(2):1456–1464, 2016.
  * [13] L. W. Liao and C. C. Yang.  On some new properties of the gamma function and the Riemann zeta function.  Math. Nachr., 257:59–66, 2003.
  * [14] L. Markus.  Differential independence of $\Gamma$ and $\zeta$.  J. Dynam. Differential Equations, 19(1):133–154, 2007.
  * [15] A. Ostrowski.  Über Dirichletsche Reihen und algebraische Differentialgleichungen.  Math. Z., 8(3-4):241–298, 1920.
  * [16] J. Steuding.  Value-distribution of $L$-functions, volume 1877 of Lecture Notes in Mathematics.  Springer, Berlin, 2007.
  * [17] S. M. Voronin.  The distribution of the nonzero values of the Riemann $\zeta$-function.  Trudy Mat. Inst. Steklov., 128:131–150, 260, 1972.  Collection of articles dedicated to Academician Ivan Matveevič Vinogradov on his eightieth birthday. II.
  * [18] S. M. Voronin.  A theorem on the “universality” of the Riemann zeta-function.  Izv. Akad. Nauk SSSR Ser. Mat., 39(3):475–486, 703, 1975.