\section{Introduction}
The square lattice Ising model is probably the most important statistical mechanical
model ever proposed in terms of its applications.

Since Onsager's \cite{Onsager} celebrated the solution of the Ising model free energy in 1944, followed by
\cite{Yang} the proof of Onsager's result for the spontaneous magnetization in 1952, physicists have devoted themselves
to the study of the problem of elucidating the susceptibility of the square lattice Ising model. As is well known, the
Onsager solution in 1944 was publicly presented without a proof \cite{ Guttmann}. Although there is no known
closed-form expression for the susceptibility $\chi$, a large body of knowledge about the susceptibility is available.
Attempting to list all these contributions is impossible in this paper.

We will, however only mention some of the research to indicate the main mathematical approaches in the investigation of
the square lattice Ising model. We refer the reader to McCoy's book \cite{McCoy} and the review article \cite{Coy} for
a review of these
developments. We crave the forgiveness of the researchers whose papers are not cited in this paper.

In 1976 Wu, McCoy, Tracy, and Barouch \cite{Wu} showed that the susceptibility could be expressed as an
infinite sum of $n$-dimensional integrals. The main result of their analysis can be given as the statement: for the
low-temperature susceptibility $(T<T_{c})$, only even powers of $n$ contribute, starting at $n = 2$, while for the
high-temperature susceptibility $(T>T_{c})$, only odd powers contribute.
In 1999, Nickel \cite{Nickel1,Nickel2} showed that the susceptibility has a natural
boundary in the complex plane and that such functions can not be $D$-finite or holonomic functions. Later it was proved
that each of these integrals $ \chi^{(n)}$ is $D$-finite but argued that their infinite sum, that is, the full
susceptibility $\chi$ itself, is not \cite{Orrick,Orrick2}. In a series of papers, Maillard and co-workers
\cite{Zenine1,Zenine2,Boukraa1,Nickel3,Bostan,Boukraa2} found the linear
ordinary differential equations (ODEs) satisfied by $\chi^{(3)},\dots, \chi^{(6)}$. In 2011, Chan and his
coworkers \cite{Chang} extended the work in Ref. \cite{Nickel1} to other two-dimensional lattices, and gave
an expansion of the scaled form of the susceptibility to unprecedented accuracy.

Almost half a century has passed since the formal expression for the magnetic susceptibility of square lattice Ising
model \cite{Coy} is obtained. There is still, however, so much discussion and uncertainty in the investigation of the
magnetic susceptibility of the square lattice Ising model.
Therefore, it is natural to ask the question: why it is so hard to obtain a closed expression for the susceptibility
in this model. To answer this question, it is proper to recall the basic definition of magnetic susceptibility, which
is given by the relation $\chi(T)=\frac{d\langle\sigma\rangle}{dh}|_{h=0}$. Here $\langle\sigma\rangle$ is the average
magnetization, $T$ denotes temperature, and $h$ is the external field. Since the average magnetization is known only in
zero magnetic fields (h=0), it is impossible to take the derivative appearing in the definition of magnetic
susceptibility.

Therefore, the susceptibility is usually studied through its relation with the zero-field spin-spin correlation function:
\begin{eqnarray*}
&&\chi(T)=\frac{1}{k_{B}NT}\sum_{i,j,k,l}\{\langle\sigma_{k,l}\sigma_{i,j}\rangle-\langle\sigma\rangle^{2}\}.
\end{eqnarray*}
Here, $i, j$ and $k, l $ denote lattice positions and run all over the system, $\langle\sigma\rangle$ is the
spontaneous magnetization. The correlation terms appearing the above general relation are approximated by replacing
the correlation function with $\langle\sigma_{0,0}\sigma_{i,j}\rangle$ or even with $\langle\sigma_{k,l}
\sigma_{i,i}\rangle$. $T_{c}$ denotes the critical temperature, and we recall that for the isotropic square
lattice Ising model, horizontal and vertical interaction constants have
the same value $J$. Then, the spontaneous magnetization is given for $T\! <\! T_{c}$ by
$\langle\sigma\rangle=(1-k_{1}^2)^{1/8}$.
Here $k_{1}=\textnormal{sinh}(2J/k_{B}T)^{-2}$ and $\langle\sigma\rangle$ is zero for $T\!>\!T_{c}$.

The correlation function appearing in the susceptibility relation
has been studied by many authors \cite{Kauf,Coy1,Coy2} and by
various methods. However, a compactly closed expression of any
$\langle\sigma_{0,0}\sigma_{i,j}\rangle$ have not been yet found
\cite{Yamada}, which indicates that the calculation of the square
lattice susceptibility is still an open problem.

Of course, applying daunting
mathematical treatments, the above-cited studies have made substantial advances in elucidating the properties of
the susceptibility of the square lattice Ising model. On the other hand, the mathematics used in these research
papers gets so complicated that it is impossible to follow up on the physical reasoning and even the mathematical
approximations used.

Therefore, it might be relevant and important to investigate the susceptibility of the square lattice Ising model
from a different perspective. To this end, we are going to exploit the recently obtained interrelation \cite{Kaya1}
for the average magnetization, which is given by the following expression as
\begin{eqnarray}
\langle\sigma\rangle=\langle\tanh[K(\sigma_{1}\cdots \sigma_{z})+H]\rangle.
\end{eqnarray}
Here, $\sigma_{1},\cdots,\sigma_{z}$ are the nearest neighbors spins around its central spin. By applying this
relationship together with a physically plausible assumption for odd spins correlation functions, we have obtained
almost exact average magnetization expressions for the 2D lattices \cite{Kaya2} and 3D cubic lattice \cite{Kaya3}
in the absence of external magnetic field.

We think it is worthwhile to mention that, in our recent papers, we have given the details of the calculations
(with some useful tricks and tools) allowing one to obtain average magnetization expressions for other 3D Ising
lattices. The obtain average magnetization expressions are almost equivalent to the exactly obtained results and
simulation data. In addition, in the calculation of the average magnetization relations, applied mathematics is
quite tractable and manageable.

Our primary goal in this paper is to calculate the susceptibility of the square lattice Ising model with the same
tractable and manageable mathematical procedures. In doing so, we would also like to gain some more insight into the
additional relevance of the conjectured assumption in our previous work. Here, the susceptibility of the square lattice
is going to be investigated by the same average magnetization interrelation and we will use the same mathematical form
of odd spins correlation function used in our previous studies. Of course, the odd spins correlation function is going
to be modified slightly to take into account the external magnetic field dependence.

To elucidate the mathematical procedure, we first calculate the susceptibility of the 1D Ising model as an example,
so that we can compare our result with the well-known exact susceptibility expression of the 1D Ising chain. In doing so,
we hope to provide adequate confirmation for the relevance of the method which will be used in the investigation of the
susceptibility of the square lattice Ising model.

This paper is organized as follows. In the next section, we are going to calculate the susceptibility of the 1D Ising
chain to elucidate and test the validity of the mathematical procedure which is going to be used in the calculation of
the susceptibility of the square Ising lattice. In the last section, the susceptibility of the square Ising lattice is
going to be calculated with the method developed for the treatment of the 1D chain. In doing so, we are going to try to
present our calculations and the necessary assumption as clearly as possible. In the same section, we are going to
discuss both the used method and some useful indications of possible future research problems which can be treated with
the method of this paper.

\section{The susceptibility of 1D Ising chain}
We start this section with the previously derived formula in reference \cite{Kaya1}, which is given by Eq. (1) in
the introduction.
To find the critical coupling strength, we need to consider $H=0$ case. So, taking into account Eq. (1) for 1D Ising
chain in the absence of external field and if the hyperbolic tangent function is expressed as
\begin{equation}
\tanh[K(\sigma_{1}+\sigma_{2})]=\frac{1}{2}(\sigma_{1}+\sigma_{2})\tanh(2K).
\end{equation}
The critical coupling strength can be readily obtained as follows. Substituting this relation into Eq. (1) leads to
$\langle \sigma\rangle= \langle\sigma\rangle \tanh(2K)$. From this equation, it is obvious that the non-zero values
for average magnetization $\langle\sigma\rangle $ is only possible at zero temperature. Here $K$ is defined as
$K=J/(k_B T)$, where $T$ is the temperature.
Thus, the critical temperature $T_{c}$ for 1D chain is equal to zero. Since the susceptibility is defined as the
derivative of the average magnetization with respect to the external magnetic field $h$, which is defined as
$H=h/(k_B T)$, we need to obtain the average magnetization of the 1D chain as a function of $H$. To this end,
we write the hyperbolic tangent function appearing in Eq. (1) as
\begin{equation}
\tanh[K(\sigma_{1}+\sigma_{2})+H]=C_{1}(\sigma_{1}+\sigma_{2})+C_{2}\sigma_{1}\sigma_{2}+C_{3}.
\end{equation}
Under the consideration of different orientations of $\sigma_{1}$ and $\sigma_{2}$, the coefficients $C_{1}$,
$C_{2}$ and $C_{3}$ can be obtained as,
\begin{eqnarray*}
&&\!\!\!\!\!\!\!\!\textstyle C_{1}(K,H)\!=\!\frac{1}{4}[\tanh(2K+H)]-\tanh(-2K+H), {}
\nonumber \\ &&
\!\!\!\!\!\!\!\!\textstyle C_{2}(K,H)\!=\!\frac{1}{4}[\tanh(2K+H)+\tanh(-2K+H)+2\tanh(H)]-\tanh(H),
{}\nonumber \\ &&
\!\!\!\!\!\!\!\!\textstyle C_{3}(K,H)\!=\!\frac{1}{4}[\tanh(2K+H)+\tanh(-2K+H)+2\tanh(H)].
\end{eqnarray*}
Substituting these final relations into Eq. (3), and then taking the average of the both side of the equation,
leads to
\begin{equation}
\langle\sigma\rangle=2C_{1}\langle\sigma\rangle+C_{2}\langle\sigma_{1}\sigma_{2}\rangle+C_{3}.
\end{equation}
Thus $\langle\sigma\rangle$ can be expressed as
\begin{equation}
\langle\sigma\rangle=\frac{C_{2}\langle\sigma_{1}\sigma_{2}\rangle+C_{3}}{1-2C_{1}}.
\end{equation}
The susceptibility is defined as
\begin{equation}
\chi(K,0)=\frac{d\langle\sigma\rangle}{dh}\Big|_{h=0} = \frac{1}{k_{B}T}\frac{d\langle\sigma\rangle}{dH}\Big|_{H=0}.
\end{equation}
After doing some algebra, $\chi(K,0)$ can be expressed as
\begin{equation}
\chi(K,0)=\frac{1}{k_{B}T}\frac{1+[\frac{1}{2}+\langle\sigma_{1}\sigma_{2}\rangle]
\textnormal{sech}^{2}(2K)}{1-\tanh(2K)}.
\end{equation}
For large values of $K$ or small values of $T$, this equation can be expressed as
\begin{equation}
\chi(K,0)=\frac{1}{k_{B}T}e^{4J/k_B T}.
\end{equation}
This final susceptibility relation is equivalent to the susceptibility relation obtained in \cite{Pathria}.
At this point it is important to point out that we have obtained the susceptibility relation without using
the well known 1D chain average magnetization relation, which is given by the following relation as
\begin{equation}
\langle\sigma\rangle=\frac{\textnormal{sinh}(H)}{[e^{-4K}+\textnormal{sinh}^{2}(H)]^{\frac{1}{2}}}.
\end{equation}

In the next section, our primary goal will be to calculate the susceptibility of the square lattice Ising model
in the same manner used in this section.
Exploiting these advantages of mathematical simplicity and tractability, we will make use of the available average
magnetization relation for the square lattice. We will work with almost the same mathematical form for the three
spins correlation function conjectured previously in reference \cite{Kaya2}. Focusing on the transition region,
we can calculate the susceptibility relation for the square lattice and we are also able to extract the values of
the magnetic critical exponent $\gamma$ for $K\!<\!K_{c}$ and $K\!>\!K_{c}$ for the square lattice.

In addition, we also give reasonable discussion for the modification of the conjectured form of the three spins
correlation function in the presence of an external magnetic field. To better understand this connection between
the three spins correlation function in the presence and the absence of an external magnetic field, we trust the
relevance of the previously obtained average magnetization relations obtained in references \cite{Kaya2,Kaya3}
in which the same mathematical form of the three spin correlation function was used.

\section{The susceptibility of the square lattice Ising Lattice}
In this section we apply the same method used in the previous section to calculate the susceptibility of the
square Ising lattice. As it is done above, the procedure starts by writing the Eq. (1) for square lattice as
\begin{equation}
\langle \sigma_{0,i}\rangle=
\langle \tanh[K(\sigma_{1,i}+\sigma_{2,i}+\sigma_{3,i} +\sigma_{4,i})+H]\rangle,
\end{equation}
where $K$ is the coupling strength and $\sigma_{0,i}$ denotes the
central spin at the $i^{th}$ site, while $\sigma_{l,i}$,
$l=1,2,3,4$, are the nearest neighbor spins around the central
spin. In what follows the index $i$ is not used since it is not necessary in the treatment. We express the
hyperbolic tangent function with the following equivalent relation as
\begin{eqnarray}
&&\tanh[K(\sigma_{1}+\sigma_{2}+\sigma_{3}+\sigma_{4})+H]=A_{1}(\sigma_{1}+\sigma_{2}+\sigma_{3}+\sigma_{4}){}
\nonumber \\ && +A_{2}(\sigma_{1}\sigma_{2}+\sigma_{2}\sigma_{3}+\sigma_{3}\sigma_{4}+\sigma_{4}\sigma_{1}+
\sigma_{1}\sigma_{3}+\sigma_{2}\sigma_{4}){}\nonumber \\ &&+A_{3}(\sigma_{1}\sigma_{2}\sigma_{3}+\sigma_{2}
\sigma_{3}\sigma_{4}+\sigma_{3}\sigma_{4}\sigma_{1}+\sigma_{4}\sigma_{1}\sigma_{2})+A_{4}
\sigma_{1}\sigma_{2}\sigma_{3}\sigma_{4}+A_{5},
\end{eqnarray}
\begin{eqnarray*}
&&\!\!\!\!\!\!\!\!\textstyle A_{1}\!=\!\frac{1}{16}[\tanh(4K\!+\!H)\!-\!\tanh(-4K\!+\!H)\!+\!2\tanh(2K\!+\!H)\!-
\!2\tanh(2K\!\!+\!\!H)],{}
\nonumber \\ &&
\!\!\!\!\!\!\!\!\textstyle A_{2}\!=\!\frac{1}{16}[\tanh(4K\!+\!H)\!+\!\tanh(-4K\!+\!H)\!-\!2\tanh(H)],{}\nonumber \\ &&
\!\!\!\!\!\!\!\!\textstyle A_{3}\!=\!\frac{1}{16}[\tanh(4K\!+\!H)\!-\!\tanh(-4K\!+\!H)\!-\!2\tanh(2K\!+\!H)\!+
\!2\tanh(2K\!\!+\!\!H)],{}
\nonumber \\ &&
\!\!\!\!\!\!\!\!\textstyle A_{4}=2A_{2}-A_{5}+\tanh(H),{}\nonumber \\ &&
\!\!\!\!\!\!\!\!\textstyle A_{5}=A_{2}+\frac{1}{4}[\tanh(2K+H)+\tanh(-2K+H)+2\tanh(H)],
\end{eqnarray*}
\begin{eqnarray}
\langle \sigma \rangle=4A_{1} \langle \sigma \rangle\!+\!A_{2}[4\langle \sigma_{1}\sigma_{2} \rangle\!+\!2\langle
\sigma_{1}\sigma_{3} \rangle]\!+\!4A_{3}\langle \sigma_{1}\sigma_{2}\sigma_{3} \rangle\!+\!A_{4}\langle \sigma_{1}
\sigma_{2}\sigma_{3}\sigma_{4} \rangle\!\!+\!\!A_{5}.
\end{eqnarray}
The three spin correlation function was conjectured as in references \cite{Kaya2,Kaya3} as
\begin{eqnarray}
\langle \sigma_{1}\sigma_{2}\sigma_{3} \rangle=a(K_{c},0)\langle \sigma \rangle+[1-a(K_{c},0)]\langle \sigma
\rangle^{\frac{1+\beta}{\beta}},
\end{eqnarray}
in the absence of external magnetic field. It is important to notice that the value of $a(K{c},0)$ was obtained
in our previous paper \cite{Kaya2}.

In this paper, we assume that the form of the three spin correlation function in the absence of external field
can also be the relevant relation for the three spin correlation function in the presence of external magnetic field. This can be expressed as,
\begin{eqnarray}
\langle \sigma_{1}\sigma_{2}\sigma_{3} \rangle=a(K,H)\langle \sigma \rangle+[1-a(K,H)]\langle \sigma
\rangle^{\frac{1+\beta}{\beta}}.
\end{eqnarray}
Thus, substituting this relation into Eq. (12) leads to
\begin{eqnarray}
f_{1}(K,H)\langle\sigma\rangle+f_{2}(K,H)\langle \sigma \rangle^{\frac{1+\beta}{\beta}}=f_{3}(K,H).
\end{eqnarray}
The functions appearing in this equation are given by
\begin{eqnarray*}
&&\!\!\!\!\!\!\!\!\textstyle f_{1}(K,H)\!=\!1-4A_{1}(K,H) -4A_{3}(K,H)a(K,H),{}
\nonumber \\ &&
\!\!\!\!\!\!\!\!\textstyle f_{2}(K,H)\!=\!-4A_{3}[1-a(K,H)],{}\nonumber \\ &&
\!\!\!\!\!\!\!\!\textstyle f_{3}(K,H)\!=\!A_{2}[4\langle \sigma_{1}\sigma_{2} \rangle+
2\langle \sigma_{1}\sigma_{3} \rangle]+A_{4}(K,H)\langle \sigma_{1}\sigma_{2}\sigma_{3}\sigma_{4} \rangle+A_{5}(K,H).
\end{eqnarray*}

Taking the derivative of the both sides of the Eq. (15) with respect to the external field $h$, the following
expression for the susceptibility of square lattice obtained as,
\begin{eqnarray}
\lim_{H\to 0} \chi(K,H)=\frac{1}{k_{B}T} \lim_{H \to 0}\frac{I_{1}(K,H)}{I_{2}(K,H)}.
\end{eqnarray}
Here $I_{1}$ and $I_{2}$ are given by the relations,
\begin{eqnarray*}
&&I_{1}(K,H)=\frac{df_{3}}{dH}-\langle \sigma \rangle \frac{df_{1}}{dH}-\langle \sigma \rangle^{\frac{1+\beta}{\beta}}
\frac{df_{2}}{dH},{}\nonumber \\&&
I_{2}(K,H)=f_{1}+\frac{1+\beta}{\beta}\langle \sigma \rangle^{\frac{1}{\beta}}f_{2},
\end{eqnarray*}
\begin{eqnarray*}
\lim_{H \to 0} \frac{df_{1}}{dH}=-\frac{1}{4}[2\tanh(4K)-4\tanh(2K)]\frac{a(K,H)}{dH}\!\!\mid_{H=0},
\end{eqnarray*}
\begin{eqnarray*}
\lim_{H \to 0} \frac{df_{2}}{dH}=\frac{1}{4}[2\tanh(4K)-4\tanh(2K)]\frac{a(K,H)}{dH}\!\!\mid_{H=0},
\end{eqnarray*}
\begin{eqnarray*}
\lim_{H \to 0} \frac{df_{3}}{dH}&=&\frac{1}{8}\tanh(4K)^{2}[4\langle \sigma_{1}\sigma_{2} \rangle \!\!\mid_{H=0}+
2\langle \sigma_{1}\sigma_{3}\rangle\!\!\mid_{H=0}] {}\nonumber\\&&+\frac{1}{8}[3-4\textnormal{ sech}(2K)^2+
\textnormal{sech}(4K)^2 ] \langle \sigma_{1}\sigma_{2}\sigma_{3}\sigma_{4} \rangle\!\!\mid_{H=0},
\end{eqnarray*}
\begin{eqnarray*}
\textstyle f_{1}(K,0)\!\!&=&1-\frac{1}{4}[4\tanh(2K)\!+\!2\tanh(4K)]\!{}\nonumber \\&&-\frac{1}{4}[2\tanh(4K)-
4\tanh(2K)]a(K,0),
\end{eqnarray*}
\begin{eqnarray*}
\textstyle f_{2}(K,0)=-\frac{1}{4}[-4\tanh(2K)+2\tanh(4K)][1-a(K,0)].
\end{eqnarray*}
For $K<K_{c}$, using the above relations, the susceptibility can be expressed as
\begin{eqnarray*}
\chi(K,0)=\frac{\frac{1}{8}\tanh(4 K)^{2}C+\frac{1}{8}[3-4 \textnormal{sech}(2K)^2+\textnormal{sech}(4K)^2 ]
\langle \sigma_{1}\sigma_{2}\sigma_{3}\sigma_{4} \rangle\!\!\mid_{H=0} }{1-[1-a(K,0)]\tanh(2K)-\frac{1}{2}[1+
a(K,0)]\tanh(4K)}.
\end{eqnarray*}
Here $C$ is defined as $C=[4\langle \sigma_{1}\sigma_{2} \rangle\!\!\mid_{H=0}+2\langle \sigma_{1}\sigma_{3}
\rangle\!\!\mid_{H=0}]$.
Now, if we can figure out the form of the unknown function $a(K,0)$, then we can calculate the susceptibility
for $K\!<\!K_{c}$.
\begin{figure*}[!hbt]
\begin{center}
\scalebox{1.1}{\includegraphics{sus1.eps}}
\end{center}
\caption{The values of $\gamma$ for the case $K<K_{c}$ in the
vicinity of $K_{c}$ for the three different values of the
arbitrary parameter $b$, which is chosen as $ b\simeq -2.68$.}
\end{figure*}
\begin{figure*}[!hbt]
\begin{center}
\scalebox{1.1}{\includegraphics{sus2.eps}}
\end{center}
\caption{The values of $\gamma$ for the case $K<K_{c}$ in the
vicinity of $K_{c}$ for the three different values of the
arbitrary parameter $b$. }
\end{figure*}
\begin{figure*}[!hbt]
\begin{center}
\scalebox{1.1}{\includegraphics{sus3.eps}}
\end{center}
\caption{The values of $\gamma$ for the case $K>K_{c}$ in the
vicinity of $K_{c}$ for the three different values of the
arbitrary parameter $b$. }
\end{figure*}
\begin{figure*}[!hbt]
\begin{center}
\scalebox{1.1}{\includegraphics{sus4.eps}}
\end{center}
\caption{The plot of $\ln(\chi)$ for $K<K_{c}$ in the vicinity of
$K_{c}$. }
\end{figure*}
\begin{figure*}[!hbt]
\begin{center}
\scalebox{1.1}{\includegraphics{sus5.eps}}
\end{center}
\caption{The plot of $\ln(\chi)$ for $K>K_{c}$ in the vicinity of
$K_{c}$. }
\end{figure*}

To this end we are going to exploit $a(K_{c},0)=0.75736$ which was
obtained in reference \cite{Kaya2}. Since our purpose is to obtain
the susceptibility function for values near to $K_{c}$, we can
expand $a(K,0)$ around $K_{c}$ as
\begin{eqnarray}
a(K,0)=a(K_{c},0)+b(K-K_{c})+\cdots,
\end{eqnarray}
here $b$ is equal to $\frac{da(K,0)}{dK}|_{K=K_{c}}$. After defining $a(K,0$), we are ready to investigate
the expression for susceptibility. The first thing one may notice is that all of the terms in the numerator
of $\chi(K,0)$ are smooth functions and they all have finite values for the any values of $K$ while the
denominators go to zero (or they assume very small values) around the critical point.

This means that the behavior of $\chi$ strongly depends on the values of denominator and it is independent of
the finite values of numerator. If the
scaling form of susceptibility, which is defined as $\chi=(1-\frac{K}{K_{c}})^{-\gamma}$, is recalled, the
critical exponent $\gamma$ can be given by the following relation as
\begin{eqnarray}
\gamma=\frac{\ln[1-[1-a(K,0)]\tanh(2K)-\frac{1}{2}[1+a(K,0)]\tanh(4K)]}{\ln(K_{c}-K)}.
\end{eqnarray}
In Fig. 1 and Fig. 2, the critical exponent $\gamma$ is plotted for different values of the parameter $b$
introduced to describe $a(K,0)$. From these figures, it is obvious that the values of $\gamma$ strongly depend
on the values of $b$ and the values of $(K-K_{c})$. It is clear from Fig. 1 that the values of $\gamma$ approach
a definite value when the values of $K$ are very close to $K_{c}$. In this limit, the value of $\gamma$ is
around $\gamma=1.7$ which is in agreement with the value $\gamma=1.75$ obtained from the scaling theory. Fig. 2 is
plotted for positive values of the parameter $b$. As seen from the figure, the values of $\gamma$ are
around $\gamma=1.1$ when $K$ is very close to $K_{c}$. Now, we are going to investigate the behavior of
the critical exponent $\gamma$ for the values of $K$ greater than $K_{c}$ in the same manner as it is treated
in the $K\!<\!K_{c}$ case. The first thing one can notice is that the terms in $I_{1}$ are smooth and they all
have finite values for all the values of $K$. This property of $I_{1}$ easily leads to the following $\gamma$ relation,
\begin{eqnarray}
\gamma=\frac{\ln[f_{1}(K,0)+9[1-\textnormal{sinh}(2K)^{-4}]f_{2}(K,0)]}{\ln(K_{c}-K)},
\end{eqnarray}
where, for the average magnetization we used the Onsager relation, $\langle\sigma\rangle=
[1-\textnormal{sinh}(2K)^{-4}]^{1/8}$, and the critical exponent $\beta$ of square lattice is equal to $1/8$.

The obtained $\gamma$ relation for the case $K\!\!>\!\!K_{c}$ is
plotted in Fig. 3. From this figure, it is obvious that the values
of $\gamma$ around the critical point $K_{c}$ do not depend on the
values of the arbitrary parameter $b$. And the values of $\gamma$
approach $0.91$ regardless of the values of the arbitrary
parameter $b$.

To elucidate the general feature of the behavior of magnetic susceptibility, the logarithm of the magnetic
susceptibility relations obtained for the cases $K\!\!<\!\!K_{c}$ and $K\!\!>\!\!K_{c}$ in this paper for the
square lattice Ising model are plotted in Fig. 4 and Fig. 5 respectively.

\section{Conclusion and Discussion}
In this paper, a new approach is introduced to investigate and calculate the susceptibility of the square
lattice Ising model. To do so, we start our investigation by applying the exact interrelation for the average
magnetization introduced previously. The application of this
interrelation in the presence of an external magnetic field produces the relations which relate average magnetization
of the lattices to the even spins and the odd spins correlation functions.

We have assumed that the even spins correlation functions behave smoothly and
continuously around the critical point. As a result, the derivative of the correlation functions with respect to
the external magnetic field gives finite values at the critical point for $H=0$. We have also assumed that the
three spin correlation function obeys the same mathematical form as conjectured in our previous works by just
replacing the constant $a$ with a function $a(K, H)$. Exploiting the previously obtained result $a(K,0)$ at
critical point, $a(K,0)$ is expanded into series around $K=K_{c}$. Since we aim to calculate the susceptibility
of the square lattice for values of $K$ near $K_{c}$, the function $a(K,0)$ can be approximated by keeping just
two terms in the expansion.

We have observed that the values of the critical exponent $\gamma$ are strongly dependent on the values of the
arbitrary parameter for the case $K<K_{c}$, but the values of $\gamma$ converge to a fixed value if the value of
arbitrary parameter takes a value around $-2.65$. Using this consideration, we have determined the values
of $\gamma=1.72$ for $K<K_{c}$. On the other hand, for the values of $K$ is greater than $K_{c}$, the values
of $\gamma$ are slightly dependent on the values of the arbitrary parameter.

In the end, $\gamma$ converges to $0.91$. We believe that the approach used in this paper is quite relevant
and important. The method can be also applicable to the more complicated untouchable problem such as calculating
the susceptibility of the 3D cubic lattice Ising model. These types of problems are going to be the subject of
our future research.


\section{Introduction}

    Measuring quantum observables is one of the bottlenecks of quantum variational algorithms.\cite{gonthier2020identifying}
    Here we consider quantum observables represented as a linear combination of Pauli products
  \begin{eqnarray}
  \HO = \sum_k c_k \HP_k,
  \end{eqnarray}
  where $c_k$ are coefficients and $\HP_k$ are tensor products of single-qubit Pauli operators.
  To measure an expectation value of $\HO$ one needs to
  transform $\HO$ into the Ising form $\HU^\dagger \HO \HU = \hat Z$, where $\HZ$ is a linear combination of Pauli products that contain only $\hat z$ operators for individual qubits. Finding such $\HU$ for an arbitrary observable is challenging
  since it is equivalent to diagonalizing $\HO$ in the computational basis. Fortunately, the expectation
  value of $\HO$ can be obtained by summing the expectation values of simpler operators that
  can be easily transformed to their Ising forms ($\HZ_n$):
  \begin{eqnarray}
  \bra{\Psi} \HO\ket{\Psi} = \sum_n \bra{\Psi} \HO_n\ket{\Psi} = \sum_n \bra{\Psi} \HU_n^\dagger \HZ_n \HU_n\ket{\Psi}.
\end{eqnarray}
   This is a partitioning approach to the measurement problem that was successfully applied in the variational quantum eigensolver (VQE),\cite{2014VQE} where $\HO$ corresponds to the system Hamiltonian, $\HH$. Even for a fixed observable
 there are usually numerous possible partitionings $\HO = \sum_n \HO_n$, where each partitioning involves its own set of
 unitary transformations $\HU_n$. For all partitionings $\HU_n$ should be straightforward to apply both on
 classical and quantum computers because one needs to know in advance $\HZ_n$ operators for classical
 post-processing and be able to apply $\HU_n$ to the quantum state $\ket{\Psi}$ before the measurement.
 The main quantity of interest that ranks partitionings in their efficiency is the total number of measurements required for achieving accuracy $\epsilon$ in estimating $\bra{\Psi} \HO \ket{\Psi}$, $N_m$. This number
 under the condition of optimal distribution of measurements between fragments $\HO_n$ is
 \begin{eqnarray}\label{eq:nm}
        N_m=\left(\frac{1}{\epsilon}\sum_n\sqrt{\text{Var}_{\Psi}(\HO_n)} \right)^{2},
 \end{eqnarray}
 where $\text{Var}_{\Psi}(\HO_n) = \bra{\Psi} \HO_n^2 \ket{\Psi} -  \bra{\Psi} \HO_n \ket{\Psi}^2$ is the quantum
 variance for each operator fragment.\cite{Crawford}

 Since one of the possible areas for demonstrating the quantum advantage is the electronic structure problem, for
 concreteness, here we focus on obtaining the expectation value of electronic Hamiltonians.
 There are multiple ways one can partition the Hamiltonian into diagonalizable fragments. For electronic
 Hamiltonians, there have been suggested two main approaches to partitioning using diagonalizable fragments in
 fermionic and qubit operator algebras. Recent comparisons of state-of-the-art techniques in the two categories shown that
 one can achieve a lower number of the measurements using techniques within the qubit operator algebra.\cite{CSA2021,Aadi_Overlapping,shlosberg2021adaptive}
 These qubit techniques are based on partitioning of the Hamiltonian to sets of commuting
 Pauli operators.\cite{yen2020measuring,MoscaA,gokhale2019minimizing,Verteletskyi:2020do,shlosberg2021adaptive}

 Due to differences in accuracy of performing one- and two-qubit transformations in
 quantum computers, it is convenient to distinguish two types of commutativity, general or full commutativity (FC) and
 a more restrictive qubit-wise commutativity (QWC). Two Pauli products qubit-wise commute only if every pair of corresponding single-qubit operators commute in them (e.g. $\hat x_1 \hat y_2 \hat z_3$ qubit-wise commutes with $\hat x_1 \hat z_3$ and does not qubit-wise commute with $\hat y_1 \hat x_2 \hat z_3$).
 Measuring a linear combination of Pauli products that qubit-wise commute require only local (one-qubit)
 Clifford transformations ($\HU_n$), while for the same task on commuting Pauli products one generally needs non-local (entangling) Clifford transformations. The main advantage of grouping schemes based on FC over QWC schemes
  is lower numbers of measurements required for estimating the Hamiltonian expectation values (\eq{eq:nm}).
  However, estimates based on \eq{eq:nm} do not account for additional uncertainties related to
  non-unit fidelities of gates that are required to introduce extra transformations needed for measurement, $\HU_n$.
  In this work, we investigate the extent to which these additional uncertainties reduce the advantage of the FC-based schemes over the QWC-based schemes in the number of measurements.

Classical shadow tomography (CST)\cite{Huang_2020} has been proposed recently as an alternative to grouping of operator terms in measurable fragments for estimating the operator expectation value.  CST techniques are based on learning about the expectation value of the operator by rotating the system with a random unitary transformations from a certain distribution and measuring the system in the computational basis. The obtained set of measurement results are projections (classical shadows) of the true system state. Expectation values of the operator for classical shadows can be calculated on a classical computer
and averaged to reproduce the expectation value of the
true system state. The Clifford group and its single-qubit subgroup distributions of unitaries have been used to draw
random unitaries.\cite{Huang_2020,Zhao_2021,Hadfield_LBCS_2006,hadfield2021adaptive,PRXQuantum:POVM}  Using the electronic structure Hamiltonians as examples of quantum observables, it was shown that for a single observable, CST schemes are generally inferior than grouping techniques in lowering the number of measurements required to estimate the expectation value.\cite{Aadi_Overlapping,shlosberg2021adaptive}
Therefore, here we consider the grouping methods to measurements but
our analysis involving Clifford transformations can be extended to the CST schemes as well.

     \section{Method}

\subsection{Error Model}

In both FC and QWC partitionings the electronic Hamiltonian in the qubit representation is
partitioned to fragment Hamiltonians, $\hat H = \sum_{n=1}^{N_f}\hat H_n$,
    \begin{equation}
         \hat H_n = \sum_{j=1}^{M_n} a_{j}P_{j},
    \end{equation}
    where $P_{j}$ are Pauli products, and $P_{j}$ within each $\hat H_n$ commute or qubit-wise commute with each other.
    For all $\hat H_n$, there are unitary transformations from the Clifford group $\hat U_n$ so that
    $\hat H_n = \hat U_n^\dagger \hat Z_n \hat U_n$.

The VQE workflow can be separated in two parts: 1) preparing the state $\ket{\Psi} = \hat U \ket{0}$
and 2) transforming $\ket{\Psi}$ to the state that will be measured: $\ket{\Phi_n} = \hat U_n \ket{\Psi}$.
The total fidelity of the $\ket{\Phi_n}$ state preparation $F_n = p q_n$ is a product of fidelities for the preparation
$\ket{\Psi}$, $p$, and of that for the additional measurement preparatory circuit $\hat U_n$, $q_n$.
FC and QWC schemes will be different in values of $q_n$ and share the $p$-part.
Since the total fidelity is multiplicative we assume that the common $p$-part
also incorporates any noise contributions appearing due to projective measurements of all qubits. For the $q_n$-part,
fidelities of one- and two-qubit gates are denoted as $f_1$ and $f_2$, respectively. Thus, if $G_{i,n}$ are the numbers of
$i$-qubit gates for the $\hat U_n$ circuit then $q_n = f_1^{G_{1,n}} f_2^{G_{2,n}}$.

    Due to non-unit gate fidelities in the state preparation we assume that instead of $\ket{\Phi_n}$
    one can only prepare the mixed state
    \begin{equation}\label{eq:rhon}
        \rho_n = F_n \ket{\Phi_n}\bra{\Phi_n} +  \frac{(1-F_n)}{d}\mathbf{1},
    \end{equation}
    where $d=2^{N_q}$ is the dimension of the $N_q$-qubit space, and $\mathbf{1}$ is the identity operator in the $d$-dimensional space.  In \eq{eq:rhon}, it is assumed that there is no systematic bias and extra gates only contribute to probability of the depolarized component.

    \subsection{Estimator characteristics}

    Here we construct an estimator for $\bra{\Psi} \hat H_n \ket{\Psi}$ assuming that we only have access to measurement
    results for the density $\rho_n$ (\eq{eq:rhon}) on the transformed fragment Hamiltonian, $\hat H_n = \hat U_n^\dagger \hat Z_n \hat U_n$
    \begin{eqnarray}
        {\rm Tr}[\hat Z_n \rho_n] &=& F_n \bra{\Phi_n} \hat Z_n \ket{\Phi_n} + (1-F_n)\text{Tr}[\hat Z_n]/d \\
        &=& F_n \bra{\Phi_n} \hat Z_n \ket{\Phi_n} = F_n \bra{\Psi} \hat H_n \ket{\Psi},
    \end{eqnarray}
    where the second equality follows from the tracelessness of $\hat Z_n$.
    Therefore, an estimator for $\bra{\Psi} \hat H_n \ket{\Psi}$, $\bar{H}_n$, can be formulated as
    \begin{eqnarray}
        \bar{H}_n = \frac{{\rm Tr}[\hat Z_n \rho_n]}{F_n}.
    \end{eqnarray}
    Its statistical variance can be obtained from the quantum variance of ${\rm Tr}[\hat Z_n \rho_n]$
    \begin{equation}
        \text{Var}(\bar{H}_n) = \frac{{\rm Tr}[\hat Z_n^2 \rho_n]-{\rm Tr}[\hat Z_n \rho_n]^2}{F_n^2}.
    \end{equation}
    Next, we would like to connect the statistical variance of our estimator with the quantum variance of the fragment
    $\text{Var}_{\Psi}(\hat H_n) = \bra{\Psi} \hat H_n^2 \ket{\Psi}- \bra{\Psi} \hat H_n \ket{\Psi}^2$.
    This connection will elucidate a functional dependence of the estimator variance on
    fidelity contributions and the operator properties. Using $\rho_n$ form in \eq{eq:rhon} gives
    \begin{eqnarray} \notag
    \text{Var}(\bar{H}_n) &=& \frac{\bra{\Psi} \hat H_n^2 \ket{\Psi}}{F_n} +\frac{1-F_n}{dF_n^2} {\rm Tr}[\hat Z_n^2] \\
    &&-\bra{\Psi} \hat H_n \ket{\Psi}^2 \\ \notag
    &=& \text{Var}_{\Psi}(\hat H_n) + \frac{1-F_n}{F_n} \bra{\Psi} \hat H_n^2 \ket{\Psi} \\
    &&+ \frac{1-F_n}{dF_n^2} {\rm Tr}[\hat Z_n^2].\label{eq:HvarP}
    \end{eqnarray}
    Both terms additional to $\text{Var}_{\Psi}(\hat H_n)$ in \eq{eq:HvarP} are positive. ${\rm Tr}[\hat Z_n^2]/d$  can be further
    simplified as ${\rm Tr}[\hat Z_n^2]/d = \sum_{j=1}^{M_n} a_j^2$ since any product of two different Pauli products is traceless.
    This gives the final expression for the estimator variance as
    \begin{eqnarray}\notag
     \text{Var}(\bar{H}_n) &=& \text{Var}_{\Psi}(\hat H_n) + \frac{1-F_n}{F_n} \bra{\Psi} \hat H_n^2 \ket{\Psi} \\
     &&+ \frac{1-F_n}{F_n^2} \sum_{j=1}^{M_n} a_j^2 .\label{eq:Hvar}
    \end{eqnarray}
 This expression allows us to formulate the minimal number of measurements for achieving $\epsilon$ accuracy with
 $66.7\%$ probability as
     \begin{equation}\label{eq:bb2}
        N_m=\left(\frac{1}{\epsilon}\sum_n\sqrt{\text{Var}(\bar{H}_n)} \right)^{2}.
    \end{equation}
Due to the relation $\text{Var}(\bar{H}_n) > \text{Var}_{\Psi}(\hat H_n)$ for $p \ne 1$ or $q_n \ne 1$, $N_m$ will always be
larger than its counterpart evaluated ignoring non-unit fidelity of gates.

To understand  why lower fidelities can be more detrimental for the QWC partitioning than for the FC one,
it is instructive to rewrite $\text{Var}(\bar{H}_n)$ as
    \begin{eqnarray}\notag
     \text{Var}(\bar{H}_n) &=& \text{Var}_{\Psi}(\hat H_n) + \frac{1-F_n}{F_n} \sum_{i\ne j} a_ia_j\bra{\Psi} \hat P_i\hat P_j \ket{\Psi} \\
     &&+(1-F_n)(F_n^{-1}+F_n^{-2}) \sum_{j=1}^{M_n} a_j^2, \label{eq:Hvar2}
    \end{eqnarray}
    where we partitioned the $\bra{\Psi} \hat H_n^2 \ket{\Psi}$ term into the diagonal ($i=j$) and off-diagonal ($i\ne j$) parts.
    For lower fidelities the last term in \eq{eq:Hvar2} will dominate because of the $F_n^{-2}$ scaling.
    In the case where $F_n$ is dominated by the uniform reduction factor $p$, the estimator variance for each
    fragment can be approximated as
         \begin{eqnarray}
     \text{Var}(\bar{H}_n) &\approx& \frac{(1-p)}{p^2} \sum_{j=1}^{M_n} a_j^2, \label{eq:HvarA}
    \end{eqnarray}
and the total number of measurements as
         \begin{eqnarray}\label{eq:peff}
         N_m \approx N_m^{(p)} = \frac{(1-p)}{\epsilon^2 p^2}\left(\sum_n\sqrt{\sum_{j=1}^{M_n} a_j^2} \right)^{2}.
    \end{eqnarray}
Clifford group transformations do not change coefficients $a_j$, hence the total sum over all $a_j^2$ in all
groups is a constant dependent on the system Hamiltonian rather than the partitioning approach.
Clearly, this constant is partitioned to smaller parts in the QWC grouping than in the FC one. It is a general
property of a square root function that $\sqrt{a+b}\le \sqrt{a}+\sqrt{b}$ for real positive $a$ and $b$. Thus, if the
total sum of $a_j^2$ is split into a larger number of smaller groups in QWC,
and QWC $N_m^{(p)}$ will be larger than that for the FC scheme.
This also follows from the fact that any QWC group is also FC group but not the other way around.

    \section{Results}
    We assess the difference in the number of measurements required in the FC and QWC schemes on a set of molecular electronic Hamiltonians. Details on the Hamiltonians are presented in Appendix \ref{sec:hamdet}. Among various heuristics for grouping commuting Pauli products\cite{Verteletskyi:2020do,yen2020measuring,tilly2021variational} we used the Sorted Insertion (SI) method\cite{Crawford} because it generally outperforms other grouping techniques in the number of measurements.\cite{CSA2021,Aadi_Overlapping}
Since measurement circuits for $\HU_i$ of the FC and QWC groups can be constructed
 using only Clifford group transformations, according to the Gottesman-Knill theorem
 these circuits can be implemented efficiently on both classical and quantum computers.\cite{gottesman1998heisenberg,10.5555/1972505}
 For the FC grouping an asymptotically optimal scaling of $O(N_q^2/\log(N_q))$ in the number of 2-qubit entangling gates and single qubit rotations can be achieved.\cite{ImprovedSimulation,patel2003efficient} We have implemented this optimal circuit synthesis and optimization method in the Tequila program\cite{2021Tequila}. The total number of CNOT, Hadamard and Phase gates for model systems averaged over measurable FC groups indeed follows
 the $O(N_q^2/\log(N_q))$ dependence according to \fig{fig:scaling}.

        \begin{figure*}[!ht]
    \center
    \includegraphics[width=1\textwidth]{Gate_Complexity.eps}
  \caption{\label{Gate Complexity} An average number of gates ($N_g$) for measurement circuits of the FC scheme as a function of $N_q^2/\log_{10}(N_q)$, where $N_q$ is the number of qubits. The averaging is done over measurable groups of the FC scheme for molecular systems described in Appendix A. The linear function $N_g = 0.72 N_q^2/\log_{10}(N_q) + 21$ (red) has $R^2 = 0.997$.}
      \label{fig:scaling}
  \end{figure*}

   Considering variety in fidelities of quantum gates available in different
    near term quantum computing architectures, we calculated $q_n$-parts of the total fidelity using two sets of 1- and 2-qubit gate fidelities. One set is based on the numbers from a recent work describing the Google superconducting quantum hardware,\cite{arute2019quantum} $f_1=0.9938$ and $f_2 = 0.9984$. This set represents what is available on a contemporary superconducting qubit device. The second set is based on a projection that in the near future we will have quantum gates with $f_1 = 0.9999$ and $f_2 = 0.999$. The $p$-component of the total fidelity depends
    not only on gate fidelities but also on the state complexity, here, we would like to see the trend by evaluating the number of
    measurement estimates for $p=1$, $p=0.8$, and $p=0.6$.

    Table \ref{tab:FCvsQWC} provides the ratios of the measurement numbers required in the FC and QWC schemes.  Ideal FC (FC-I) has always lower $N_m$ than that of QWC (QWC-I) except for H$_2$. The ratio between the numbers of measurements  grows with the system size and reaches more than 20 for the largest system (N$_2$).
    When we introduce uncertainty related to non-unit fidelity of gates, the ratios decrease but in most cases indicate
    superiority of the FC scheme. This reduction in ratios is the consequence of generally lower fidelities
    ($F_n$) in FC groups compared to those of QWC groups because of a larger number of gates needed for
    measurements in the former. For small systems, H$_2$ and LiH, non-unit fidelity of gates lead to lower numbers of measurements for the QWC scheme (ratios $<1$). Interestingly, accounting for the reduction of the circuit fidelity
    due to the $\ket{\Psi}$ state preparation ($p=0.8$ and $p=0.6$) increases the ratios for the small systems and even
    reverses the trend for H$_2$ and for LiH with the second set of gate fidelities. This detrimental effect of $p$ reduction on
    QWC $N_m$'s is attributed to the contribution of \eq{eq:peff} that favours a smaller number of groups. The effect from
    \eq{eq:peff} generally competes with other terms in \eq{eq:Hvar2} and is clearly seen in cases
    of H$_2$, LiH, and NH$_3$ (gate fidelities of a contemporary device). For the other systems, $p$ reduction leads to reduction
    of the ratios.

       \begin{table*}[!htbp]
        \setlength\tabcolsep{0pt}
        \caption{Ratios of the number of measurements $N_m$ to achieve $\epsilon=1 \times 10^{-3}$ Hartree accuracy in the energy expectation value (with 66.7\% probability) assuming the optimal distribution of measurements for individual fragments.
        Two partitioning methods are considered: FC and QWC. Three gate fidelity sets are used: ideal, unit fidelities (-I); fidelities available on a contemporary device,   set 1 (-C); and near future fidelities, set 2 (-F). }
            \centering
        {\begin{tabular*}{1\textwidth}{@{\extracolsep{\fill}}  l l l l l l l l}
        \toprule
         Systems & H$_2$ & LiH & BeH$_2$ & H$_2$O & NH$_3$ & N$_2$ \\
        $N$ & 4 & 12 & 14 & 14 & 16 & 20 \\
        \midrule
     \multicolumn{7}{c}{$p=1$} \\
          $N_m$(QWC-C)/ $N_m$(FC-C) & 0.965 & 0.602 & 1.91 & 1.15 & 1.20 & 1.67 \\
         $N_m$(QWC-F)/ $N_m$(FC-F)  & 0.993 & 0.990 & 4.54 & 2.72 & 4.17 & 14.1 \\
         $N_m$(QWC-I)/ $N_m$(FC-I) & 1.00 & 1.08 & 5.51 & 3.24 & 5.30 & 23.2 \\
        \midrule
        \multicolumn{7}{c}{$p=0.8$} \\
         $N_m$(QWC-C)/ $N_m$(FC-C) & 1.04 & 0.912 & 1.37 & 1.14 & 1.18 & 1.24 \\
         $N_m$(QWC-F)/ $N_m$(FC-F) & 1.06 & 1.11 & 1.68 & 1.36 & 1.75 & 2.12 \\ \midrule
          \multicolumn{7}{c}{$p=0.6$} \\
         $N_m$(QWC-C)/ $N_m$(FC-C) & 1.09 & 0.979 & 1.29 & 1.15 & 1.20 & 1.20 \\
         $N_m$(QWC-F)/ $N_m$(FC-F) & 1.10 & 1.16 & 1.55 & 1.33 & 1.65 & 1.85 \\
        \bottomrule

        \end{tabular*}
        }
        \label{tab:FCvsQWC}
    \end{table*}

    \section{Conclusions}

Using a simple unbiased depolarization error model we have derived an estimation
for the number of measurements needed in the VQE measurement schemes involving local and non-local qubit transformations.
Considering realistic fidelities for current and near future superconducting quantum processors and using the implementation of the non-local unitary transformations minimizing the number of CNOT gates, it was found that the non-local scheme based on grouping commuting Pauli products generally requires  fewer measurements than the corresponding local scheme based on grouping of qubit-wise commuting Pauli products.
Due to a larger number of additional gates the advantage of the non-local scheme decreases with reduction in gate fidelities. Another trend is that the advantage from the FC scheme is diminished even further if one accounts for the fidelity reduction
due to the state preparation circuit. However, lowering of the state preparation fidelity can also exhibit an interesting opposite trend that favors the FC scheme because it has a smaller number of groups.
Overall, we found that for all systems larger than LiH (in the STO-3G basis) the number of measurements
needed in the FC scheme is lower than that in the QWC scheme even with accounting for
all uncertainties introduced by non-unit gate fidelities.

    \section*{Acknowledgements}
     A.F.I. acknowledges financial support from the Google Quantum Research Program, Early Researcher Award, and Zapata Computing Inc. This research was enabled in part by support provided by Compute Ontario and Compute Canada.


\section{Introduction}

Spoken language identification (LID) is the task of automatically identifying the language of an utterance. Speech-based LID is used as a pre-processing step in several applications, such as automatic call routing, multilingual spoken translation and human-machine communication systems, multilingual speech transcription systems and spoken document retrieval. SLR is also often used in the area of intelligence and security.

In order to encourage the research on multilingual phenomena and advance the development of language and dialect
recognition technologies, the Oriental Language Recognition (OLR) Challenge has been organized annually since 2016~\cite{olr2016,olr2017,olr2018,olr2019,olr2020}. The sixth OLR challenge, denoted by OLR 2021 Challenge \cite{wang2021olr}, included two LID tasks and two ASR tasks, both evaluating either constrained or unconstrained training scenarios.  The constrained  LID task is a cross-domain closed-set identification task  with 13 target languages. Only the data provided by the organizers can be used to build the system. The unconstrained LID task is a closed-set identification task with 17 languages. Here,  utterances obtained from real-life environments are used for evaluation.  Any data and models are allowed for system training and development. As in NIST LRE15, the OLR 2021 Challenge uses $C_{\rm avg}$ as the principle evaluation metric.

The Tallinn University of Technology (TalTech) team participated in both language identification tasks. We relied on transfer learning: in the constrained task, we first trained a multilingual automatic speech recognition (ASR) model and finetuned it for language recognition, similarly as proposed in \cite{wang21z_interspeech}. In the unconstrained task, we finetuned the multilingual XLSR-53 wav2vec2.0 model \cite{conneau2020unsupervised} first on the VoxLingua107 dataset \cite{valk2021slt} and then on the target language/dialect data.

\section{Previous work}

Various pretraining and transfer learning schemes have been a popular approach to LID for a long time. In the phonotactic approach, one first trains a multilingual \cite{matejka2005phonotactic} or several monolingual phone recognizers \cite{zissman1996comparison} on transcribed training data. Streams of phone symbols are then recognized from test utterances. Phone stream from a test utterance is finally classified using language-specific phone n-gram models.

Another successful pretraining strategy is based on bottleneck features (BNFs) \cite{matejka2014neural,jiang2014deep,fer2015multilingual}.  Transcribed monolingual or multilingual spoken data is used for training a DNN-based acoustic model that has a narrow bottleneck layer close to the output layer. After training, the layers after the bottleneck layer are discarded and the output from the bottleneck layer is used as a frame-level feature extractor. Bottleneck features provide impressive performance improvements with both i-vector as well as neural network based frontends \cite{snyder2018spoken}.  The best-scoring submission to the OLR 2020 Challenge tasks 1 and 3 was also based on bottleneck features, extracted from an end-to-end multilingual Conformer-based ASR model \cite{duroselle2021language}.

Several works have shown that finetuning a pretrained feature extractor with the rest of the backend model for spoken language identification gives additional benefits with regard to using a frozen model for extracting bottleneck features. In \cite{wang2021end}, the shared encoder of a multilingual Conformer-based ASR model was finetuned for dialect identification, achieving strong results in the second task of the OLR 2020 Challenge. Recently, self-supervised models trained on large amounts of untranscribed data have been reported to achieve large improvements on language recognition. In \cite{tjandra2021improved} it was shown that multilingual wav2vec2.0 model XLSR-53 can be finetuned for language recognition using very little training data. XLS-R \cite{babu2021xlsr}, an improved version of the XLSR-53 wav2vec2.0 model, achieved state-of-the-art results on the VoxLingua107 development set.

\section{Task 1: Constrained Language Identification}

Task 1 of the OLR 2021 Challenge is a closed-set identification task with 13 target languages.  Training and test data is recorded in various different  environments, making it a cross-domain identification task. Only speech data provided by the organizers is allowed to be used for training. Non-speech data from external sources is allowed for data augmentation purposes.

\subsection{Data}

\begin{table}[tb]
\caption{Task 1 training and development data.}
\label{tab:task1-data}
\begin{tabular}{l|rr|rr|r}
\toprule
         & \multicolumn{2}{l|}{All training data} & \multicolumn{2}{l|}{Transcr. data} & Dev \\
\midrule
Language         & \#Utts             & \#h            & \#Utts            & \#h          & \#Utts        \\
\midrule
Kazakh    & 10748            & 18.2           & 4192            & 9.5           & 0           \\
Hokkien   & 9305             & 16.5           & 5163            & 15.6          & 1998        \\
Shanghainese & 9472             & 12.8           & 5374            & 12.1          & 1796        \\
Sichuanese  & 9487             & 11.5           & 5378            & 10.9          & 1798        \\
Tibetan    & 19480            & 18.0           & 5378            & 11.9          & 0           \\
Uyghur    & 13497            & 25.0           & 5382            & 13.6          & 0           \\
Cantonese    & 18572            & 21.5           & 5705            & 7.6           & 4353        \\
Indonesian    & 16053            & 18.9           & 5733            & 7.5           & 1800        \\
Japanese    & 17728            & 15.4           & 5734            & 5.8           & 4116        \\
Korean    & 17597            & 16.8           & 8800            & 6.0           & 3725        \\
Russian    & 14581            & 20.6           & 8967            & 9.9           & 3744        \\
Vietnamese    & 16902            & 21.0           & 8982            & 8.4           & 1800        \\
Mandarin    & 19302            & 22.0           & 11094           & 7.6           & 3600       \\
\bottomrule
\end{tabular}
\end{table}

We used all the provided training data for the 13 target languages for training the models, except for the OLR 2020 test data, that was used as a development set. A subset of the training data comes with transcripts. This data was used for pretraining the ASR model. The amount of data for each target language is summarized in Table~\ref{tab:task1-data}.
The background noises in the Freesound portion of the MUSAN corpus \cite{musan2015} and  simulated small, medium and large room impulse responses \cite{ko2017study} were used for data augmentation.

\subsection{Models}

In the constrained task, our system was a combination of four models. One of the models was a Resnet-style model, trained on the provided training from scratch. The other three were finetuned from Conformer-based multilingual ASR models, with different pooling methods and training criteria.

\subsubsection{Resnet-style model}

\begin{figure}[tbh]
\centering
\includegraphics[width=0.6\linewidth]{Augmix1.pdf}
\caption{{\it Sample augmentation pipeline in AugMix (largely based on \cite{hendrycks2019augmix}). The number and type of transforms in each augmentation path is randomly sampled. The augmented signals are mixed using randomly sampled $w_1$ and $w_2$. The mixed augmented signal is finally interpolated with the original signal, using a sampled weight $m$.}}
\label{fig:augmix1}
\end{figure}

The Resnet-style model is derived from the x-vector paradigm~\cite{snyder2018x,snyder2018spoken}, with several enhancements. During training, we apply on-the fly data augmentation using AugMix \cite{hendrycks2019augmix}, by randomly distorting the training data using a mix of reverberation and noise augmentation. A clean training segment is cloned into several copies. A different randomly drawn series of augmentation transforms, each possibly with random parameters, is applied serially to each of the copies. Then, the augmented copies are mixed with each other (using randomly sampled weights) and the resulting super-augmentated sample is finally mixed with the original clean sample, using a randomly sampled interpolation coefficient (see Figure \ref{fig:augmix1}).  The benefit of this method, compared to using pre-generated static augmentations, is that there is a lot of variety in the training data: each training sample is a result of several random transforms, applied in random order and with a random weight.

For frame-level feature extraction, we use the Resnet34 \cite{cai2018exploring, he2016deep} architecture where the basic convolutional blocks with residual connections are replaced with squeeze-and-excitation modules \cite{hu2018squeeze,zhou2019deep}.
The statistics pooling layer that maps frame-level features to segment level features is replaced in our model with a multi-head attention layer \cite{bahdanau2014neural} that has been shown to provide superior performance  \cite{zhu2018self, rahman2018attention, okabe2018attentive, safari2019self}. From among many variants of multi-head attention used in previous studies, we employ the one described in \cite{zhu2018self}: frame level representations are first mapped to $N_{\rm att}$ outputs ($N_{\rm att}=128$ in our model), using a $1\times1$ convolution and a ReLU nonlinearity; from this representation, each attention head (we used  $N_{\rm heads}=5$ heads) computes it's own softmax-based weight distribution over the input utterance; finally, weighted mean and standard deviation are computed over the frame level features for each head, resulting in $N_{\rm heads} \times 512 \times 2$ segment-level representations.

The structure of the embedding model is summarized in Table \ref{tab:architecture}. The variables $F$ and $T$ refer to the number of filterbanks and the number of time frames in the utterance. In all experiments, we used $F=30$.

\begin{table}[tb]
\caption{The neural network architecture of the Resnet model. Notation: \textit{SE/res} -- squeeze-and-attention block and residual connections; \textit{F} -- number of filterbanks in the input features; \textit{T} -- number of frames of input; \textit{FC} -- fully connected layer. }
\label{tab:architecture}
\footnotesize
\begin{tabular}{lllc}
\toprule
\textit{Layer}          &\textit{Spatial Size}        & \textit{\#Channels} &   \textit{Kernel}       \\
\midrule
Input          & $F \times T$            & 1              & -            \\
\multicolumn{4}{l}{\it{Frame level representations}} \\
Pre-resnet     & $F \times T$       & 64             & 7 $\times$ 7 \\
Res-block 1  & $F/2 \times T/2$   & 64       & $3 \times \begin{bmatrix} 3 \times 3 \\ 3 \times 3 \\ SE/res \end{bmatrix}$  \\
Res-block 2 & $F/4 \times T/4$   & 128        & $4 \times \begin{bmatrix} 3 \times 3 \\ 3 \times 3 \\ SE/res \end{bmatrix}$ \\
Res-block 3 & $F/8 \times T/8$   & 256       & $6 \times \begin{bmatrix} 3 \times 3 \\ 3 \times 3 \\ SE/res \end{bmatrix}$ \\
Res-block 4 & $F/16 \times T/16$ & 512        & $3 \times \begin{bmatrix} 3 \times 3 \\ 3 \times 3 \\ SE/res \end{bmatrix}$ \\
Post-resnet & $1 \times T/16$    & 512            & $F/16  \times 1$ \\
\multicolumn{4}{l}{\it{Segment-level representations}} \\
Pooling   & 1    & $5 \times 512 \times 2$ &  Attent. stats \\
Embedding        & 1    & 512     &  Dense \\
FC        & 1    & 512         &   Dense \\
Output       & 1    & \#Languages         &   Softmax \\
\bottomrule
\end{tabular}
\end{table}

The models are implemented in PyTorch \cite{NEURIPS2019_9015} using a framework developed in our lab.

\subsubsection{Conformer ASR model}

The Conformer-based ASR model was trained on pooled provided training data that came with transcripts. As a development set, we used the OLR 2020 test data. We applied a number of text normalization steps to the transcripts before training: for Cantonese, Mandarin, Mandarin dialects and Japanese, all whitespace symbols were deleted. For all languages, all punctuation symbols were deleted. For Kazakh, the Arabic script was transliterated into Cyrillic.

The ASR model uses a byte-pair encoding (BPE) vocabulary of 20\,000 units, shared over all languages. This is much larger than a typical BPE vocabulary size for a monolingual model. However, since the BPE tokens are shared over 13 languages, it averages to around 1500 BPE token per language which is typical for low resource end-to-end ASR models. The multilingual ASR model is an encoder-decoder based model that uses Conformer as an encoder and Transformer as a decoder. Some important hyperparameters of the model are listen in Table \ref{tab:conformer}. The model was trained on speed-perturbed data and SpecAugment \cite{park2019specaugment} (without warping) was applied during training. The number of training epochs was 40. The model was trained on four GPUs, with a dynamic batch size of 30M frames per GPU, which translates to 4620 updates per epoch.
After every epoch, the model's performance on development data was measured, and the final model was averaged over 10 best-performing models.

\begin{table}[tb]
\caption{Hyperparameters of the Conformer ASR model}
\label{tab:conformer}
\centering
\begin{tabular}{lr}
\toprule
\textit{Conformer encoder} &  \\
\midrule
Number of blocks & 12 \\
Linear dimensionality & 2048 \\
Dropout rate & 0.1 \\
Output size & 256 \\
Num. attention heads & 8 \\
\midrule
\textit{Transformer decoder} & \\
\midrule
Linear dimensionality & 2048 \\
Number of blocks & 6 \\
Num. attention heads & 4 \\
\midrule
\textit{Training} & \\
\midrule
CTC weight & 0.3 \\
Label smoothing & 0.1 \\
\bottomrule
\end{tabular}
\end{table}

After training of the ASR model, the encoder part of this model was taken as the backbone of the language recognition model. The encoder's outputs were fed through a pooling layer. We experimented with different pooling methods: attentive statistics pooling \cite{okabe2018attentive}, multi-head attention (MHA) pooling  \cite{india2019self}, global multi-head attention (GMHA) pooling \cite{wang2020multi}. Two fully connected layers with the ReLu non-linearity and BatchNorm were appended to the pooling layer. As a training criterion, we experimented with cross-entropy (CE) loss and additive angular margin (AAM) loss \cite{deng2019arcface}. The ASR encoder part of the model was trained together with the rest of the model, without any learning rate scaling. For some models, stochastic weight averaging (SWA) \cite{izmailov2018averaging} was applied during the last 30\% of training.

The language recognition model was trained on all available provided training data for the 13 target languages. Language embeddings were extracted from the first fully-connected layer after the pooling layer. The dimensionality of the embeddings was 512.

The Conformer-based model for ASR was trained using ESPnet \cite{watanabe2018espnet}. For pooling operations and backend scoring (see below), ASV Subtools \cite{tong2021asv} was used.

\subsubsection{Back-end modeling}

For backend scoring, we used a multinomial regression model. Language embeddings were normalized and mean-centered based on the training data. No LDA was used in this task. The logistic regression model was rebalanced using priors proportional to the inverse amount of data for each language, in order to remove any language bias from the model.

We combine the scores of various systems using calibrated combination weights. For finding the model combination weights, we optimize  the parameters of a linear model based on the log-likelihood cost metric (CLLR) on the development data, using L-BFGS as the optimizer. Our own Pytorch-based calibration implementation was used which is freely available\footnote{\url{https://github.com/alumae/sv_score_calibration}}.

\subsection{Results}

Our results for Task 1, together with a comparison to the official baseline and the winning team's results, are listed in Table~\ref{tab:task1-res}. Note that the Progress set is a 30\% subset of the Test set, hence the results should be very similar. However, some teams improved their system after submission to Progress set leaderboard was closed which explains the large difference.
It can be seen that transfer learning from an ASR model gives a massive boost to the system performance. Our best single model is based on the trained ASR Conformer encoder, used global multihead attention pooling, cross-entropy loss and stochastic weight averaging. Fusion of four models gives only a slight improvement both on development and progress sets. The fusion was also used in our official test set submission.

\begin{table*}[tb]
\centering
\caption{Results on Task 1 with various systems and their combination.}
\label{tab:task1-res}
\begin{tabular}{lll|ll|ll|ll}
\toprule
              &                 &           & \multicolumn{2}{c|}{\textit{Dev}} & \multicolumn{2}{c|}{\textit{Progress}} &   \multicolumn{2}{c}{\textit{Test}} \\
\textit{Backbone}      & \textit{Pooling}         & \textit{Criterion} & \textit{$C_{\rm avg}$} & \textit{EER}  & \textit{$C_{\rm avg}$}  & \textit{EER} & \textit{$C_{\rm avg}$}           & \textit{EER}          \\
\midrule
Resnet        & Attentive stats & CE        & 0.0515      & 6.90     & 0.0601        & 5.72        \\
ASR Conformer & MHA             & CE        & 0.0110       & 1.68     & 0.0101        & 1.15         \\
ASR Conformer & GMHA            & CE + SWA  & 0.0081      & 1.27     & 0.0080         & 0.91       \\
ASR Conformer & MHA             & AAM       & 0.0129      & 2.19     & 0.0133        & 1.68        \\
\midrule
\multicolumn{3}{l|}{\textbf{Fusion}}                  & 0.0078      & 1.20     & \textbf{0.0074}        &\textbf{0.86}  & 0.0079 & 0.86    \\
\midrule
Baseline \cite{wang2021olr} &  & &  & & 0.0826 & 9.04  & 0.0817 & 8.98  \\
Team X-Voice (ranked 1st) \cite{lyuant2021}  &  & &  & & \textbf{0.0062} &	\textbf{0.65} & 0.0025 & 0.27 \\
Team funspeech (ranked 3rd) \cite{wangfunspeech2021} & & &  & & 0.0080 & 0.91 & 0.0083 & 0.93 \\
\bottomrule
\end{tabular}
\end{table*}

\section{Task 2: Unconstrained Language Identification}

Task 2 of the challenge is also a closed-set language identification task. However, in this task the number of target languages is 17 and the test data is more challenging, originating from real-life environments. In this task, any data, including externally available pretrained models, are allowed to be used for system building.

\subsection{Data}

\begin{table}[tbh]
\centering
\caption{Task 2 training and development data. For languages marked with an asterisk, the data originates from Mozilla CommonVoice.}
\label{tab:task2-data}
\begin{tabular}{l|rr|r}
\toprule
         & \multicolumn{2}{l|}{Training data} & Dev  \\
\midrule
Language & \#Utts             & \#Hours       & \#Utts        \\
\midrule
Hokkien (Minnan)   & 9305             & 16.5           & 1998 \\
Shanghainese & 9472             & 12.8           & 1796 \\
Sichuanese  & 9487             & 11.5           & 1798 \\
Tibetan    & 19480            & 18.0           & 0    \\
Uyghur    & 13497            & 25.0           & 0    \\
Cantonese    & 18572            & 21.5           & 4353    \\
English$^\star$       & 15000            & 20.5           & 3000 \\
Hindi$^\star$       & 6812             & 8.1            & 2020 \\
Indonesian    & 16053            & 18.9           & 1800 \\
Japanese    & 17728            & 15.4           & 4116 \\
Korean    & 17597            & 16.8           & 3725 \\
Malay$^\star$       & 15000            & 26.2           & 1800 \\
Russian    & 14581            & 20.6           & 3744 \\
Telugu$^\star$       & 4004             & 5.1            & 444  \\
Thai$^\star$       & 15000            & 18.2           & 3000 \\
Vietnamese    & 16902            & 21.0           & 1800 \\
Mandarin    & 19302            & 22.0           & 3600 \\
\bottomrule
\end{tabular}
\end{table}

For training the backend and for final finetuning of the language embedding model, we constructed a custom training  set covering all 17 languages (see Table \ref{tab:task2-data}). For most languages, we used the data provided in the OLR 2021 Challenge training dataset. However, for English, Hindi, Malay, Telugu and Thai, there is no training data provided. For those languages, we used data from Mozilla CommonVoice \cite{ardila2019common}.  For languages with a lot of CommonVoice data, we limited the data to a random 15000 utterance subset.  The development set was compiled from OLR2020 test data and CommonVoice data for the languages not covered by the OLR data. We acknowledge that the development set does not reflect the characteristics of the test data very well, since it mostly contains read speech, while the test data originates from diverse real-life situations.

For training the language embedding models, we relied a lot on the recently released VoxLingua107 dataset. The VoxLingua107 dataset \cite{valk2021slt}\footnote{Available at \url{http://bark.phon.ioc.ee/voxlingua107/}} is a large-scale dataset for training spoken language identification models that work well on diverse real-life data.

VoxLingua107 was compiled from automatically scraped YouTube data.
The data collection process consists of the following steps: first, semi-random trigram search phrases were generated from the Wikipedia text corpus of the particular language. The search phrases were used to retrieve YouTube videos whose title or description matched the search phrase. Text-based language  identification was used for filtering out the videos with the title and description likely not in the given language. Audio tracks of the videos were downsampled to 16\,kHz.
Speech activity detection and speaker diarization were applied for extracting segments from the videos that contain speech. Long speech segments were split into utterance-like subsegments of up to 20 seconds in length.
Data-driven post-filtering was used to remove segments from the database that were likely not in the given language, increasing the proportion of correctly labeled segments in the dataset to 98\%, based on crowd-sourced verification. Some numerical facts about the VoxLingua107 training data are given in Table \ref{tab:voxlingua107}.

Most of the Task 2 target languages are covered by VoxLingua107. The exceptions are the Mandarin dialects, Tibetan, Uyghur and Cantonese.

\begin{table}[tb]
\caption{Statistics about the VoxLingua107 dataset.}
\label{tab:voxlingua107}
\centering
\begin{tabular}{lr}
\toprule
Number of languages  & 107 \\
Total number of videos & 64110 \\
Total number of hours & 6682 \\
Average number of hours per language & 62 \\
Average number of utterances per language & 23709 \\
Total amount of audio (uncompressed, in GB) & 758 \\
\bottomrule
\end{tabular}
\end{table}

\subsection{Models}

\begin{table*}[tb]
\caption{Results on Task 2 with various systems.}
\label{tab:task2-res}
\centering
\begin{tabular}{l|ll|ll|ll}
\toprule
              & \multicolumn{2}{c|}{\textit{Dev}} & \multicolumn{2}{c|}{\textit{Progress}} & \multicolumn{2}{c}{\textit{Test}}\\
\midrule
\textit{Model}    & \textit{$C_{\rm avg}$}       & \textit{EER}        & \textit{$C_{\rm avg}$}        & \textit{EER}         & \textit{$C_{\rm avg}$}        & \textit{EER}    \\
\midrule
Resnet, trained on VoxLingua107                                  & 0.053      & 5.15      & 0.078         & 7.98        \\
Resnet, trained on VoxLingua107, finetuned on training data      & 0.015      & 1.53      & 0.055         & 5.27        \\
XLSR-53, finetuned on VoxLingua107                        & 0.017      & 1.76      & 0.016         & 1.69        \\
XLSR-53, finetuned on training data                       & 0.003      & 0.29      & 0.044         & 3.78        \\
\textbf{XLSR-53, finetuned on VoxLingua107, then on training data} & 0.006      & 0.78      & 0.012         & 0.93   & 0.012 & 0.94     \\
\midrule
Team X-Voice (ranked 1st)  \cite{lyuant2021}  & & & 0.017 & 1.85 & 0.004 & 0.42 \\
Team nisp\_speech (ranked 3rd) \cite{jinnisp2021}  & & &        &      & 0.032 & 3.23 \\
\bottomrule
\end{tabular}
\end{table*}

We experimented with using the XLSR-53 wav2vec2.0 model~\cite{conneau2020unsupervised} as the backbone of our language embedding model. XLSR-53 is a large pretrained model trained on unlabeled multilingual data. The model is trained by jointly solving a contrastive task over masked latent speech representations  and learning a quantization of the latents shared across languages.  The model contains a convolutional feature encoder that maps raw audio to latent speech representations which are fed to a Transformer network that outputs context representations. XLSR-53 is pretrained on
56\,000 hours of speech data from 53 languages.

We used XLSR-53 as follows: the outputs from the wav2vec2 model were fed through an attentive pooling  layer, a fully connected layer with ReLU and BatchNorm, and the final output layer, corresponding to the languages of the training set. During training, the learning rate corresponding to the XLSR-53 model was set to 0.01 times lower than the base learning rate. We experimented with three finetuning scenarios: using VoxLingua107, using the training data of the 17 languages (OLR + CommonVoice), and finetuning first on the VoxLingua107 data, followed by the training data for the 17 languages.

We also experimented with a Resnet model, trained on VoxLingua107, and then finetuned on the OLR training data.

The same on-the-fly data augmentation strategy as was used in Task 1 was also used here for all models.

\subsection{Back-end modeling}

As in Task 1, we used a multinomial regression model for backend scoring. Language embeddings were normalized and mean-centered based on the training data. The resulting embeddings were reduced to 50-dimensional vectors using LDA. No rebalancing of the model bias was used in this task, since it was found to slightly hurt the system performance on the progress set.

\subsection{Results}

The results of various systems are listed in Table \ref{tab:task2-res}. It can be seen that using XLSR-53 as the basis for finetuning models gives a big boost to the system performance, particularly on the  progress set. Our internally compiled development set had quite different performance trends, compared to the progress set. This is probably due to the fact that the development set contained clean dictated speech, whereas the progress set contained data ``from the wild''.

The best performing-model on the progress set was the XLSR-53 model, finetuned first on VoxLingua107 and then on the training data that covered all 17 target languages and dialects. This system was also used in our official test set submission.

\section{Post-evaluation analysis}

In this section, we analyze various aspects of pretraining for spoken language identification.

Since at the time of writing this paper, the OLR 2021 Challenge test data for Task 2  was not available, we performed the experiments on Task 1 and Task 2 of OLR 2020. Note that rules of the OLR 2020 Challenge restrict the participants to use only the provided training data that includes data from the previous OLR challenges and the THCHS30 Chinese speech corpus \cite{wang2015thchs} which thereby also rules out using transfer learning. Therefore, we don't claim any state-of-the art results on OLR 2020 test data as we deviate from the official rules.

Task 1 of the OLR  2020 Challenge is a closed-set identification task with six target languages (Cantonese, Korean, Indonesian, Japanese, Vietnamese, Russian). The utterances are recorded via different recording equipments and environments, making it a cross-channel LID task. Similarly to the official baseline model, we used the provided training data of the OLR  2020 Challenge (without the dialect data) for language enrollment. We omitted the datasets that were used as development test data in Task 1 or Task 2. In Task 1, the development test data was the Task 2 data of AP19-OLR (AP19-OLR-channel). This is consistent with the official baseline system developed for the OLR 2020 Challenge.

Task 3 of OLR 2020 Challenge is a closed-set identification task under noisy conditions. The language of each utterance is one of the five target languages (Cantonese, Korean, Japanese, Mandarin Chinese and Russian). The enrollment data for this task is all of the OLR 2020 Challenge training data for the particular five languages.

As the third and a very different experimental dataset, we used the Lwazi corpus \cite{feld2009multilingual}. It contains audio samples from the 11 ofﬁcial languages in South Africa: Afrikaans, English, isiNdebele, isiXhosa, isiZulu, Sepedi, Sesotho, Setswana, Siswati, Tshivenda, and Xitsonga \cite{feld2009multilingual}. This data
was collected in South Africa over the telephone and is stored using 8 kHz sample rate. Many of the utterances in the corpus are very short which makes the data challenging for LID models.
The test and training data for the experiments was composed as follows: first, we selected 150 random speakers from the dataset, and included all the utterances from those speakers into the test set. The training set was composed by randomly selecting 2 hours of speech per language while excluding the utterances of the speakers in the test set. This training and test split is similar as was done in \cite{feng2019low}, except that we use a larger test set. For transfer learning experiments, the data was resampled to 16 kHz.

\subsection{What to use as a pretrained model?}

\begin{table*}[]
\caption{Comparison of different pretrained models on three datasets. In all cases, the resulting model is used without finetuning on the target training data.}
\label{tab:pretraining}
\begin{tabular}{l|rr|rr|rr}
\toprule
                                                                          & \multicolumn{2}{l|}{\textit{OLR 2020 Task 1}} & \multicolumn{2}{l|}{\textit{OLR 2020 Task 3}} & \multicolumn{2}{c}{\textit{Lwazi}}        \\
\midrule
\textit{Embeddings}        & \textit{$C_{\rm avg}$}              & \textit{EER}             & \textit{$C_{\rm avg}$}              & \textit{EER}             & \textit{$C_{\rm avg}$}   & EER   \\
\midrule
Official baseline (no pretraining) \cite{olr2020, feng2019low}                                       & 0.1321            & 14.58           & 0.0715            & 7.14            &        & 35.30 \\
Our baseline (Resnet + MHA, no pretraining)                                               & 0.1105            & 11.87           & 0.0709            & 6.80            & 0.2572 & 26.36 \\
Resnet, trained using VoxLingua107                                        & 0.0249            & 2.76            & 0.0147            & 1.33            & 0.1185 & 12.41 \\
WenetSpeech ASR Conformer enc., finetuned on VoxLingua107              & 0.0071            & 0.83            & 0.0043            & 0.39            & 0.0947 & 9.68  \\
CV multiling. ASR Conformer enc., finetuned on VoxLingua107 & 0.0061            & 0.59            & 0.0089            & 0.78            & 0.0985 & 10.13 \\
XLSR-53 wav2vec2.0, finetuned to VoxLingua107                                        & 0.0052            & 0.56            & 0.0057            & 0.57            & 0.1319 & 13.73 \\
XLS-R-300M  wav2vec2.0, finetuned on VoxLingua107                                     & \textbf{0.0021}            & \textbf{0.21}            & \textbf{0.0029}            & \textbf{0.28}            & \textbf{0.0809} & \textbf{8.30} \\
\bottomrule
\end{tabular}
\end{table*}

In this section, we compare different pretrained models as starting points for training a language embedding model on VoxLingua107. We experimented with the following models:
\begin{itemize}
    \item Resnet model, randomly initialized (i.e., no pretraining before training with VoxLingua107)
    \item The encoder part of the Conformer-based sequence-to-sequence ASR model\footnote{\url{https://huggingface.co/espnet/pengcheng_guo_wenetspeech_asr_train_asr_raw_zh_char}} trained on WenetSpeech \cite{zhang2021wenetspeech}, a 10K+ hour corpus of Mandarin, containing speech from YouTube and podcasts.
    \item The encoder part of the Conformer-based  sequence-to-sequence multilingual ASR model\footnote{\url{https://huggingface.co/espnet/Shinji_Watanabe_open_li52_asr_train_asr_raw_bpe7000_valid.acc.ave}} trained on the Mozilla CommonVoice corpus \cite{ardila2019common}, covering 52 languages.
    \item XLSR-53, the multilingual wav2vec2.0 model \cite{tjandra2021improved}.
    \item XLS-R-300M, an improved version of XLSR-53. It is trained on half a million hours of publicly available speech data in 128 languages, including the VoxLingua107 corpus.
\end{itemize}
Results on the three test datasets are listed in Table \ref{tab:pretraining}. In all cases, the resulting model is used for extracting embeddings for the language enrollment data and the final classifier is built using LDA and logistic regression.
We did not finetune the models on the target language enrollment data, except in the ``no pretraining'' case where only the enrollment data was used for training the embedding model. As we already experienced with the OLR 2020 Challenge Task 3 results, transfer learning from ASR encoder or wav2vec2.0 models results in large gains, compared to the Resnet model trained on VoxLingua107 from scratch. It is surprising to see that a monolingual Conformer model trained on WenetSpeech results in comparable results to the multilingual ASR model and XLSR-53. This suggests that having diverse real-world data among the pretraining data is as important as multilinguality, as WenetSpeech contains speech ``from the wild'' while CommonVoice has only prompted speech. Finally, we see that the new multilingual wav2vec2.0 model gives the best results across all datasets.  We also see that pretraining is much more effective on the OLR data than on the Lwazi corpus. This can be probably explained by two factors: first, only Afrikaans and English among the 11 target languages in Lwazi are covered by VoxLingua107. Second, all the pretrained models are trained on (mostly) 16 kHz data, making them less effective for 8 kHz telephone speech.

\subsection{VoxLingua107: does size or variety matter?}

\begin{figure}[tb]
  \centering
  \includegraphics[width=0.90\linewidth]{ablation.pdf}
  \caption{Language identification performance when the amount of pretraining data is reduced, by either removing random utterances (but keeping all the languages) or by removing whole languages.}
  \label{fig:ablation}
\end{figure}

Using the full VoxLingua107 dataset for pretraining resulted in substantial improvements throughout the experiments. Here we analyze the effect of the size and language variety of the pretraining data on transfer learning performance. This is done by training Resnet-based language embedding models on various subsets of VoxLingua107. In order to factor out external effects, we do not use pretrained Conformer and wav2vec2.0 models here.

We created several subsets of the VoxLingua107 data, each containing either 50\%, 25\% or 12.5\% of the original data. The subsetting was performed in two different ways: either per utterance or per language. In the first case, a certain proportion of utterances were kept while keeping the number of languages intact, while in the second case, all utterances of randomly selected languages were kept, so that the number of utterances would be roughly  the same as in the first case.

Figure \ref{fig:ablation} shows the EER on OLR 2020 Task 3 and Lwazi data, when using the described subsets of VoxLingua107 for pretraining. No finetuning was performed in those experiments. It is clear from the experiments that the number of different languages present in the pretraining data is more important than the number of utterances.

In the case of OLR Task 3 experiments, the fast performance drop when pretraining languages are removed can be partly explained by the fact that the language subsets drop some of the OLR Task 3 target languages (Cantonese, Korean, Japanese, Mandarin Chinese and Russian). Indeed, the 50\% language subset includes only Korean, and the even smaller subsets contain none of the target languages. Therefore, we trained another model on only the four target languages (except Cantonese) of VoxLingua107 and used this as the embedding extractor. This resulted in EER of 8.1\% which is worse than the model trained on the 14 language (12.5\%) subset of VoxLingua107 where none of the target languages is present (EER of 5.7\%). This suggests that although having the target languages among the pretraining languages is beneficial, it's more important to pretrain on a large variety of languages.

\subsection{How much data do we need for training the backend?}

We investigated how much enrollment data is needed for training the logistic regression based backend model with the best embedding model based on XLS-R-300M and VoxLingua107. For this, we created random subsets of OLR 2020 Task 1 and Task 3 training data with increasing number of utterances per target language. Figure \ref{fig:utts_vs_eer} shows that with just one utterance per language, this approach already outperforms the baseline model that doesn't use pretraining, but uses around 10\,000 training utterances per language. With only five utterances per language, we already achieve EER rates of below 1\%.

\begin{figure}[tb]
  \centering
  \includegraphics[width=\linewidth]{utts_vs_eer.pdf}
  \caption{Number of enrollment utterances per language \textit{vs} EER. Note the log axes.}
  \label{fig:utts_vs_eer}
\end{figure}

\section{Conclusion}

We described the TalTech spoken language identification systems for the OLR 2021 Challenge.  For the constrained task, our best-performing single system was trained via transfer learning from the encoder part of a Conformer-Transformer based multilingual model. For the unconstrained task, we relied on two important external resources: the XLSR-53 pretrained multilingual wav2vec2.0 model, and the VoxLingua107 corpus. Our best performing model was trained via transfer learning from XLSR-53 using two finetuning steps: first using VoxLingua107, followed by target language training data.

Post-evaluation experiments showed that both ASR encoder models as well as self-supervised wav2vec2.0 models serve as good starting points for training a language embedding extractor on the VoxLingua107 dataset.
We also showed that the variety of languages in VoxLingua107 is more important than the total amount speech. Finally, we demonstrated that on cross-channel and noisy test sets of OLR 2020, using pretrained embeddings requires only five language enrollment utterances to achieve below 1\% EER.

\bibliographystyle{IEEEbib}