\section*{Introduction}
Ever since 1994 when Seiberg and Witten  
derived exact low-energy Wilsonian effective action of (pure) $SU(2)$ ${\cal N}=2$ SYM \cite{Seiberg:1994rs},
 the interest in this kind of theories has been remaining extremely high. The reason is their remarkably rich physical and 
mathematical content. In fact, these theories provide a framework to address in a precise manner such 
problems as strong coupling, non-perturbative effects and confinement in non-Abelian   
gauge theory (so relevant, for instance, in the Standard Model). The impact of Seiberg-Witten theory 
in pure mathematics is also very substantial. Likely, the most famous applications are in algebraic geometry and topology of four-dimensional 
differentiable manifolds where e.g. the notions of Seiberg-Witten and Gromov-Witten invariants are 
of primary importance. 

The effective action of ${\cal N}=2$ SYM is given in terms of prepotential: a holomorphic 
function of the vacuum expectation values (VEV) of the vector multiplet scalar field. 
Large VEV expansion of the prepotential reveals its structure as sum of classical, one-loop 
and instanton contributions. Many researchers tried to restore the instanton contributions 
directly from the microscopic theory, but they succeeded only in the case of the first few instantons \cite{Dorey:2002ik}. Actual progress has been achieved with the idea of using equivariant 
localization techniques in the moduli space of instantons \cite{Flume:2001nc}, 
\cite{Flume:2001kb}, 
especially in combination 
with the introduction of the so-called $\Omega$ background (see \cite{Nekrasov:2002qd} and 
further developments \cite{Flume:2002az}, \cite{Nekrasov:2003rj}, \cite{Bruzzo:2002xf}). Considering theory in $\Omega$-background effectively embeds the system 
in a finite volume $\sim \frac{1}{\epsilon_1\epsilon_2}$, where the parameters 
$\epsilon_1$, $\epsilon_2$ are sort of 
angular velocities on orthogonal planes of (Euclidean) 4d space-time. This makes 
the partition function a finite, well defined quantity (commonly referred as Nekrasov 
partition function). Then the corresponding free energy coincides with (generalized) 
prepotential. The usual SW prepotential is recovered simply by sending the parameters 
$\epsilon_{1,2}\rightarrow 0$. Another, crucial consequence of introducing this background 
is the fact that instanton moduli integrals are localized at finitely many points. 
This property eventually leads to an elegant combinatorial formula for instanton 
contributions \cite{Flume:2002az}.

Later developments are even more surprising. It appears that introduction of $\Omega$ 
background is not merely a regularization trick. Thus,  keeping $\epsilon_{1,2}$ finite 
a deep relation between conformal blocks of 2d CFT and Nekrasov partition function 
\cite{Alday:2009aq}
emerges, so that the Virasoro central charge is related to these parameters, the masses of 
hypermultiplets specify inserted primary fields, while VEVs identify 
the states of the intermediate channel.

The special case  $\epsilon_{1}=-\epsilon_{2}$ bridges the theory with topological string, 
$\epsilon$-expansion of Nekrasov partition function coinciding with topological (string) genus expansion. 

Another special case of great interest is the Nekrasov-Shatashvili (NS) limit \cite{Nekrasov:2009rc}
when one of 
parameters, say $\epsilon_1$ is kept finite while $\epsilon_2=0$. From AGT point of 
view this case corresponds to semiclassical CFT when the central charge 
$c\rightarrow \infty$. Besides this, another interesting link to quantum integrable 
system emerges, now the remaining nonzero parameter $\epsilon_1$ being related to 
Plank constant. In NS limit many quantities familiar from original Seiberg-Witten 
theory become deformed or quantised in rather simple manner. 

In particular the algebraic equation defining 
Seiberg-Witten curve, becomes a finite difference equation \cite{Poghossian:2010pn}, 
which in terms of 
related integrable system is nothing but Baxter's $TQ$ equation (for a later development see also 
\cite{Bourgine:2017aqk}).
 Through discrete Fourier 
transform one gets a linear differential equation \cite{Fucito:2011pn}, which from 2d CFT 
perspective is 
the null vector decoupling equation \cite{Belavin:1984vu} in the semiclassical limit. 
This relation was an object of intensive investigations in the last decade (see e.g.    
\cite{Alday:2009fs,Mironov:2009uv,Maruyoshi:2010iu,Marshakov:2010fx,Nekrasov:2013xda,
Piatek:2011tp,Ashok:2015gfa,Poghossian:2016rzb,Poghosyan:2016mkh,
Nekrasov:2017gzb,Jeong:2017mfh}).
 

More recently, moving from Gaiotto's idea of looking at these
 equations as quantum versions of the (suitable power of the) SW differential \cite{Gaiotto:2009we}\footnote{In other words, the SW differential gives way to the oper upon quantisation.}, it has been proposed to investigate their monodromies (quantum periods over cycles) through the connection (Stokes) multipliers appearing in the  ODE/IM correspondence \cite{Dorey:1998pt,Bazhanov:1998wj}: \cite{Fioravanti:2019vxi} describes the general idea by exemplifying it in the simple case of pure $SU(2)$ gauge theory\footnote{Very few details are given for the cases with matter in the fundamental.} and in particular the link between the $a$-period and the Baxter transfer matrix $T$ function. In this perspective, Thermodynamic Bethe Ansatz  (TBA)-like considerations about pure $SU(2)$ gauge theory were initiated in \cite{Gaiotto:2014bza} at zero modulus (of the Coulomb branch) for the dual period $a_D$, and then more recently pursued by \cite{Grassi:2019coc}.



In fact, in this paper we show how to compute the gauge A-periods of the pure $SU(3)$ theory (without any matter hypermultiplet: {\it cf.} \cite{Klemm:1994qs} and \cite{Argyres:1994xh} for what concerns the generalisation of SW theory to higher rank gauge groups) as Floquet monodromy coefficients of the aforementioned differential equation (in the complex domain). Then, we propose a connexion between them and the integrable Baxter's $T$ function which extends non trivially what happens in the $SU(2)$ case and shows that the latter is not an accident. More in details, we obtain a third order linear differential equation with some similarities (and differences\footnote{The main relevant difference is, as in simpler $SU(2)$ case \cite{Zamolodchikov:2000unpb} and \cite{Fioravanti:2019vxi}, the exit of the oper parameter ($M>1/2$ in \cite{Dorey:1998pt}) from the range of validity with the appearance of a extra irregular singularity in zero (besides that at $\infty$).}) with the third order oper of  ODE/IM correspondence in  
\cite{Dorey:1999pv}, \cite{Bazhanov:2001xm}. As the latter correspond to some 'minimal' case $M>1/2$, we may conjecture, along the lines of \cite{Dorey:1999uk} and \cite{Zamolodchikov:2000unpb}, that we are describing $A_2$-Toda CFT with central charge $c=98$.


In a very interesting unfinished paper \cite{Zamolodchikov:2000unpb} Alexei Zamolodchikov has proposed ODE/IM for the Liouville CFT TBA. Special attention has been payed to the self-dual case $c=25$, when the related ODE becomes the {\it modified} Mathieu equation and a elegant 
relationship between Floquet exponent and Baxter's $T$ function has been suggested. As written, the implication of this conjecture for the period on the $A$-cycle of (effective) $SU(2)$ gauge theory has been highlighted and used by \cite{Fioravanti:2019vxi}. But it was not clear from there if and how it is possible to generalise this beautiful connection between transfer matrix and periods for higher rank groups.  

In the details of this paper we derive $QQ$ and $TQ$ functional relations (see eqs. (\ref{Qc_QQ}), (\ref{T_Q})) and extend Zamolodchikov's 
conjecture for the case of gauge group $SU(3)$ (see eq. (\ref{conj_su3})). We show that numerical integration of the differential equation leads to 
evaluation of the 'quantum' $A$-cycle periods $a_{1,2,3}$ at any point of the moduli 
space of vacua parametrised by the vector multiplet scalar VEV's 
$u_2=\langle \textbf{tr}\,\phi^2\rangle$ and $u_3=\langle \textbf{tr}\,\phi^3\rangle$ 
even for large values of $q$ at which the instanton series diverges.    
We have checked that the numerical results at small $q$ are in excellent agreement with instanton 
calculation. Thus the main message of this paper is that the differential equation provides 
an excellent tool for investigation of deformed SW theory in its entire range 
from weak to strong coupling.    

The paper is organized as follows:

Section \ref{NPF} is a short review on instanton calculus for $SU(N)$ SW theory without
hypers in $\Omega$-background. 
Here one can find explicit expressions as a sum over (multiple) Young diagrams for 
Nekrasov partition function and VEV's  $\langle \textbf{tr}\,\phi^J\rangle$.
 
Section \ref{BDE} is a brief introduction to deformed SW theory. 
We present the main results of \cite{Poghossian:2010pn}  in a form convenient 
for our present purposes.  
Starting from section \ref{DE} we consider the case of $SU(3)$ theory. The main 
tool of our investigation, a third order linear ODE is derived and its asymptotic 
solutions are found.

In section \ref{FR} we identify a unique solution $\chi(x)$ which rapidly vanishes for 
large negative values of the argument $x\rightarrow -\infty $. The three quantities 
$Q_{1,2,2}$ are defined as coefficients of expansion of $\chi(x)$ in terms of 
three independent solutions $U_{1,2,3}(x)$ defined in asymptotic region $x\gg 0$. 
Investigating symmetries of the differential equation we find a system of difference 
equations for $Q_k$ and their analogs $\bar{Q}_k$ obtained by flipping the sign of
parameter $u_3\rightarrow -u_3$. Based on this $QQ$ system we introduce Baxter's 
$T$ function and write down corresponding $TQ$ relations. 
    
In section \ref{a_cycles} we show how numerical integration of the 
differential equation along imaginary direction with standard boundary 
conditions allows one to find the monodromy matrix and corresponding 
Floquet exponents, which in the context of gauge theory, coincide with the $A$-cycle 
periods $a_{1,2,3}$. We have convincingly demonstrated the correctness of this identities
 trough comparison with instanton computation. But the main value of this method is 
 that it makes accessible also the region of large coupling constants, which is beyond 
 the reach of instanton calculus. Eventually, we close this section by suggesting a simple relation between Baxter's $T$-function and 
$A$-cycle periods $a_{1,2,3}$ of $SU(3)$ theory, which can be thought of as a natural extension of Alexei Zamolodchikov's 
conjecture relating Floquet exponent of Mathieu equation to Baxter's $T$ 
function in $c=25$ Liouville CFT.

Finally appendix \ref{TQ_proof} contains few technical details 
for derivation of the $TQ$ relation.  

           
  

\section{Nekrasov partition function and the VEVs $\langle \textbf{tr}\,\phi^J\rangle $}
\label{NPF}
Consider pure $SU(N)$ theory without hypers in $\Omega$-background.
The instanton part of partition function is given by \cite{Nekrasov:2002qd}
\begin{eqnarray} 
Z_{inst}(\mathbf{a},\epsilon_1,\epsilon_2,q)=\sum_{\vec{Y}}Z_{\vec{Y}} \left((-)^Nq\right)^{|\vec{Y}|}
\,,
\label{z_inst}
\end{eqnarray}
where sum runs over all $N$-tuples of Young diagrams $\vec{Y}=(Y_1,\cdots,Y_N)$ 
, $|\vec{Y}|$ is the total number all boxes,
$\mathbf{a}=(a_1,a_2,\cdots,a_N)$ are VEV's of adjoint scalar from 
${\cal N}=2$ vector multiplet, $\epsilon_1$, $\epsilon_2$, as already mentioned, 
parametrize the $\Omega$-background  and the instanton counting parameter $q=\exp 2\pi i 
\tau $, with  $\tau=\frac{i}{g^2}+\frac{\theta}{2\pi}$
being the (complexified) coupling constant. The coefficients $Z_{\vec{Y}}$ are factorized 
as
\begin{eqnarray} 
\label{ZY_factorization}
Z_{\vec{Y}}=\prod_{u,v=1}^N\frac{1}{P(Y_u,a_u|Y_v,a_v)}\,\,,
\end{eqnarray} 
where the factors $P(\lambda,a|\mu,b)$ for 
arbitrary pair of Young diagrams $\lambda,\mu$ and associated VEV parameters $a$, $b$, 
are given explicitly by the formula \cite{Flume:2002az} 
\begin{eqnarray}
\label{P_factor} 
&&P(\lambda ,a|\mu ,b)=\\
&&\quad \prod_{s\in \lambda}(a-b+\epsilon_1(1+L_\mu(s))-\epsilon_2 A_\lambda(s))
\prod_{s\in \mu}(a-b-\epsilon_1L_\lambda(s)+(1+\epsilon_2 A_\lambda(s)))\nonumber
\end{eqnarray}
If one specifies location of a box $s$ by its horizontal and vertical 
coordinates $(i,j)$, so that $(1,1)$ corresponds to the corner box, its
leg length $L_\lambda(s)$ and arm length $A_\lambda(s)$ with respect to the 
diagram $\lambda$ ($s$ does not necessarily belong to $\lambda$) are defined as
\begin{eqnarray} 
A_\lambda(s)=\lambda_i-j; \qquad\qquad L_\lambda(s)=\lambda_j'-i\,,
\end{eqnarray}
where $\lambda_i$ ($\lambda_j'$) is $i$-th column ($j$-th row) of diagram $\lambda$ 
with convention that when $i$ exceeds the number of columns 
($j$ exceeds the number of rows) of $\lambda$, one simply sets $\lambda_i=0$
($\lambda_j'=0$).
The instanton part of (deformed) prepotential is given by \cite{Nekrasov:2002qd}
\begin{eqnarray} 
F_{inst}({\bf a},q)=-\epsilon_1\epsilon_2\log Z_{inst}\,.
\end{eqnarray}
Instanton calculus allows one to obtain also the VEV's $\langle \textbf{tr}\,\phi^J\rangle$, 
$\phi$ being the adjoint scalar of vector multiplet:
\begin{eqnarray}
\langle \textbf{tr}\,\phi^J\rangle=\sum_{i=1}^Na_u^J+
Z_{inst}^{-1}\sum_{\vec{Y}}Z_{\vec{Y}}{\mathcal O}_{\vec{Y}}^J q^{|\vec{Y}|}
\,,
\label{tr_phiJ}
\end{eqnarray} 
where $Z_{\vec{Y}}$ is already defined by (\ref{ZY_factorization}), (\ref{P_factor}), 
and \cite{Losev:2003py,Flume:2004rp}
\begin{eqnarray} 
{\mathcal O}_{\vec{Y}}^J=\sum_{u=1}^N\,\sum_{(i,j)\in Y_u}
\left(\left(a_u+\epsilon_1 i+\epsilon_2(j-1)\right)^J+
\left(a_u+\epsilon_1 (i-1)+\epsilon_2j\right)^J\right.\nonumber\\
-\left.\left(a_u+\epsilon_1 (i-1)+\epsilon_2(j-1)\right)^J
-\left(a_u+\epsilon_1 i+\epsilon_2 j\right)^J\right). 
\end{eqnarray}   
\section{A Baxter difference equation}
\label{BDE}
\subsection{Bethe ansatz equation for NS limit}
It was shown in \cite{Poghossian:2010pn} that in NS limit $\epsilon_2\rightarrow 0$, the sum 
(\ref{z_inst}) is dominated by a single term corresponding to a unique 
array of Young diagrams $\vec{Y}^{(cr)}$ specified by properties 
(the $i$-th column length of a diagram $Y_u$ will be denoted as $Y_{u,i}$):
\begin{itemize}
\item{Though the total number of boxes $\rightarrow \infty $ in 
$\epsilon_2\rightarrow 0$ limit the rescaled column lengths 
$\epsilon_2 Y^{(cr)}_{u,i}$, converge to finite values 
\[
\xi_{u,i}=\lim_{\epsilon_2\rightarrow 0}\epsilon_2 Y^{(cr)}_{u,i}\,.
\]
}
\item{
The rescaled column lengths  at small $q$ behave as $\xi_{u,i}\sim O(q^i)$. 
This means in particular, that in order to achieve accuracy up to $q^L$, 
it is consistent to consider restricted Young diagrams with number of columns $\le L$. 
}
\item{
Up to arbitrary order $q^L$ the quantities 
\[
x_{u,i}=a_u+\epsilon_1(i-1)+\xi_{u,i}
\]
satisfy the Bethe-ansatz equations (for each $u=1,2,\cdots N$)
\begin{eqnarray}
-q \prod_{v,j}^{N,L} \frac{(x_{u,i}-x_{v,j}-\epsilon_1)(x_{u,i}-x_{v,j}^0+ 
	\epsilon_1)}{(x_{u,i}-x_{v,j}+\epsilon_1)(x_{u,i}-x_{v,j}^0-\epsilon_1)}
=\prod_{v=1}^N(x_{u,i}-a_v+\epsilon_1)(a_v-x_{u,i})\,,\quad
\label{BA}
\end{eqnarray}
where, by definition
\[
x_{u,i}^0=a_u+\epsilon_1(i-1)\,.
\]
}	
\end{itemize}
The system of equations (\ref{BA}) together with the property $\xi_{u,i}\sim O(q^i)$ 
uniquely fixes the quantities $x_{u,i}$ up to order $q^L$. Of course, calculations become 
more cumbersome if one increases $L$. Examples of explicit computations for first few 
values of $L$ can be found in \cite{Poghossian:2010pn}. 
\subsection{Baxter's difference equation and deformed Seiberg-Witten 'curve'}
The BA equations can be transformed into a difference equation  \cite{Poghossian:2010pn}
\begin{eqnarray}
Y(z+\epsilon_1)+\frac{q}{\epsilon_1^{2N}}Y(z-\epsilon_1)=
\epsilon_1^{-N} P_N(z+\epsilon_1)\, Y(z)\,,
\label{difference_eq}
\end{eqnarray} 
where $Y(z)$ is an entire function with zeros located at $z=x_{u,i}$:
\begin{eqnarray}
Y(z)=\prod_{u=1}^N e^{\frac{z}{\epsilon_1} \psi(\frac{a_u}{\epsilon_1})}\prod_{i=1}^\infty 
\left(1-\frac{z}{x_{u,i}}\right)e^{z/x_{u,i}^0}\,,
\label{Y}
\end{eqnarray}
and 
\[
\psi(x)=\frac{d}{dx}\log \Gamma(x)
\]
is the logarithmic derivative of Gauss' gamma-function. Finally $P_N(z)$ is 
an $N$-th order polynomial which parametrises the Coulomb branch of the theory. 
Explicit expressions of coefficients of this polynomial in terms of VEVs 
\begin{eqnarray}
u_J\equiv \langle {\textbf tr} \phi^J \rangle
\label{u_J_def}
\end{eqnarray}  
will be presented later for the case of our current interest $N=3$. For more 
general cases one can refer to \cite{Poghossian:2010pn}. Now, let us briefly recall how the 
difference equation (\ref{difference_eq}) is related to the Seiberg-Witten curve.
Introducing the function 
\[
y(z)=\epsilon_1^{N}\,\,\,
\frac{Y(z)}{Y(z-\epsilon_1)}
\]
one can rewrite (\ref{difference_eq}) as
\begin{eqnarray}
y(z)+\frac{q}{y(z-\epsilon_1)}=P_N(z)\,.
\label{DSW} 
\end{eqnarray}
At large $z$ the function $y(z)$ behaves as
\[
y(z)=z^N(1+O(1/z))\,.
\] 
Notice that setting $\epsilon_1=0$ in (\ref{DSW}) one obtains 
an equation of hyperelliptic curve, which is just the Seiberg-Witten curve. 
When $\epsilon_1\neq 0$, everything goes surprisingly similar to the original Seiberg-Witten theory. For example the r\^ole of Seiberg-Witten 
differential is played anew by the quantity 
\[
\lambda_{SW}=z \frac{d}{dz}\log y(z) \, ,
\]  
and, as in the undeformed theory, the expectation values (\ref{u_J_def}) are given by the contour integrals
\begin{eqnarray}
\langle \textbf{tr}\,\phi^J\rangle =\oint_{\cal C} \frac{dz}{2\pi i} z^J\partial_z \log y(z)\,,
\label{u_J_int}
\end{eqnarray}
where $\cal C$ is a large contour, enclosing all zeros and poles of $y(z)$.
\subsection{Details on $SU(3)$ theory}
Without any essential loss of generality, from now on we will assume that 
\begin{eqnarray}
u_1 \equiv  \langle\textbf{tr}\,\phi \rangle=a_1+a_2+a_3=0\,.
\end{eqnarray}
Representing $y(z)$ as a power series in $1/z$
\begin{eqnarray}
y(z)=z^3(1+c_1z^{-1}+c_2z^{-2}+c_3z^{-3}+\cdots )
\label{y_expansion}
\end{eqnarray} 
and inserting in eq. (\ref{u_J_int}) one easily finds the relations 
\begin{eqnarray} 
c_1=0; \qquad c_2=-\frac{u_2}{2}; \qquad c_3=-\frac{u_3}{3}\,.
\label{c_u}
\end{eqnarray} 
Now, consistency of (\ref{y_expansion}), (\ref{c_u}) and (\ref{DSW}) 
immediately specifies the polynomial $P_3(z)$ (we omit the 
subscript $3$, since only the case $N=3$ will be considered later on)
\begin{eqnarray} 
P(z)=z^3-\frac{u_2}{2}\,z-\frac{u_3}{3}\,\,.
\end{eqnarray}
\section{The differential equation and its asymptotic solutions}
\label{DE}
\subsection{Derivation of the differential equation}
To keep expressions simple, from now on we will set $\epsilon_1=1$. 
In fact, at any stage the $\epsilon_1$ dependence can be easily restored on dimensional grounds. Taking the results of previous subsection, the difference equation
for $N=3$ case (\ref{difference_eq}) can be rewritten as
\begin{eqnarray}
Y(z)-\left(z^3-\frac{u_2}{2}z-\frac{u_3}{3}\right) Y(z-1)+q\,Y(z-2)=0
\,,
\label{difference_eq_3}
\end{eqnarray}
By means of inverse Fourier transform, following \cite{Fucito:2011pn,
Nekrasov:2013xda,Poghossian:2016rzb}, from (\ref{difference_eq_3}) we can derive a third order linear differential equation for the function
\begin{eqnarray} 
f(x)=\sum_{z\in\mathbb{Z}+a}e^{x (z+1)}Y(z)\,.
\label{f_series}
\end{eqnarray}
At least when  $|q|$ is sufficiently small, it is expected that the series 
is convergent for finite $x$, 
provided $a$ takes one of the three possible values $a_1$, $a_2$ or $a_3$.  
Taking into account the difference relation (\ref{difference_eq_3}), one 
can easily check that  the function (\ref{f_series}) solves the 
differential equation
\begin{eqnarray} 
-f^{\,'''}(x)+\frac{u_2}{2}\,f^{\,'}(x)+\left(e^{-x}+q\,e^x+\frac{u_3}{3}\right)
f(x)=0\,.
\label{diff_eq1}
\end{eqnarray} 
Denoting
\[ 
q=\Lambda^6 
\]
and shifting the variable
\[
x\rightarrow x-\log \Lambda^3
\] 
the differential equation (\ref{diff_eq1}) may be cast into a more symmetric 
form
\begin{eqnarray} 
-f^{\,'''}(x)+\frac{u_2}{2}\,f^{\,'}(x)+\left(\Lambda^3(e^x+e^{-x})+\frac{u_3}{3}\right)
f(x)=0\,.
\label{diff_eq2}
\end{eqnarray}
\subsection{Solutions at $x\rightarrow \pm \infty $}
Physics leads us to introduce parameters $p_1$, $p_2$, $p_3$ 
satisfying $p_1+p_2+p_3=0$ such that
\begin{eqnarray}
u_2=p_1^2+p_2^2+p_3^2=2(p_1^2+p_2^2+p_1p_2); \quad 
u_3=p_1^3+p_2^3+p_3^3=-3p_1 p_2 (p_1+p_2)\,, \qquad
\label{u_2_3}
\end{eqnarray}
as in the weak coupling limit $\Lambda \rightarrow 0$ the parameters $p_i$ and $a_i$, respectively, coincide. 

At large positive values $x\gg 3\ln \Lambda$ the term $\Lambda^3 e^{-x}$ in (\ref{diff_eq2}) can be neglected. In this region the differential equation 
can be solved in terms of hypergeometric function $_0F_2(a,b;z)$ 
defined by the power series
\begin{eqnarray} 
_0F_2(a,b;z)=\sum_{k=0}^{\infty}\frac{z^k}{k!(a)_k(b)_k}\,,
\end{eqnarray}
where 
\begin{eqnarray}
(x)_k=x(x+1)\cdots(x+k-1)
\end{eqnarray}
is the Pochhammer symbol. Three linearly independent solutions can be chosen as
\begin{eqnarray}
U_i(x)\approx e^{(x+3\theta )p_i}\,\,
_0F_2(1+p_i-p_j,1+p_i-p_k;e^{x+3\theta })\,\,,
\label{U_solutions}
\end{eqnarray} 
where by definition 
\begin{eqnarray}
\Lambda \equiv \exp \theta
\end{eqnarray}
and the indices $(i,j,k)$ 
are cyclic permutations of $(1,2,3)$. We used the symbol $\approx$ 
in (\ref{U_solutions}) to mean that the approximations of the solutions hold, striktly speaking, only for $x\gg 3\theta $ (at leading order). In the end, we must verify that the Wronskian of the three solutions (\ref{U_solutions}) (below and later on, for brevity, we use the notation $p_{ij}\equiv p_i-p_j$)
\begin{eqnarray}
Wr[U_1(x),U_2(x),U_3(x)]\equiv \det \left( 
\begin{array}{ccc}
U_1(x)&U_2(x)&U_3(x)\\
U_1^{'}(x)&U_2^{'}(x)&U_3^{'}(x)\\
U_1^{''}(x)&U_2^{''}(x)&U_3^{''}(x)
\end{array}
\right) 
=p_{12}p_{23}p_{31}\quad 
\label{U_wr}
\end{eqnarray}
is not zero provided the parameters $p_i$ are pairwise different. Thus, (\ref{U_wr}) confirms that generically the $U_i(x)$ are linearly independent 
and constitute a basis in the space of all solutions. 

Similarly in region $x\ll -3\theta $ the term $\Lambda^3 e^{x}$ 
of (\ref{diff_eq2}) becomes negligible and one can write down the three linear independent solutions 
\begin{eqnarray}
V_i(x)\approx e^{(x-3\theta )p_i}\,\,
_0F_2(1-p_i+p_j,1-p_i+p_k;-e^{-x+3\theta })\,\,.
\label{V_solutions}
\end{eqnarray}
In fact, we obtain the same result for the Wronskian \begin{eqnarray}
Wr[V_1(x),V_2(x),V_3(x)] =p_{12}p_{23}p_{31} \, .
\label{V_wr}
\end{eqnarray}
\section{The functional relations}
\label{FR}
\subsection{The $QQ$ relations}
All three solutions $V_i(x)$ grow very fast at $x\rightarrow -\infty$, but 
there is a special linear combination (unique, up to a common constant factor) which vanishes in this limit. If it is the fastest one (as we suspect), this solution is usually referred to as subdominant. Using formulae for asymptotics of $_0F_2$, which 
can be found e.g. in \cite{NIST:DLMF}, we are able to establish that the correct combination is
 \begin{eqnarray}
\label{chi}
&\chi(x)=
\frac{ \Gamma (p_{12}) \Gamma ( p_{13}) }{4 \pi ^2}\, V_1(x)+
\frac{ \Gamma (p_{23}) \Gamma ( p_{21}) }{4 \pi ^2}\, V_2(x)+
\frac{ \Gamma (p_{31}) \Gamma ( p_{32}) }{4 \pi ^2}\, V_3(x)\,.
\end{eqnarray}
Its asymptotic expansion at $x\rightarrow -\infty$ is given by
 \begin{eqnarray}
\chi (x)=
\frac{v^{-\frac{1}{3}} e^{-3 v^{1/3}} }{2 \pi  \sqrt{3}}
\left(
1-\left(\frac{1}{9}-\frac{u_2}{2}\right) v^{-\frac{1}{3}}+ 
\left(\frac{u_2^2}{8}-\frac{5 u_2}{36}+\frac{u_3}{6}+\frac{2}{81}\right)v^{-\frac{2 }{3}} \right.\nonumber\\-
\left. \left(-\frac{u_2^3}{48}+\frac{u_2^2}{18}-\frac{u_3 u_2}{12}-\frac{13 u_2}{324}+\frac{7 u_3}{54}+\frac{14}{2187}\right)v^{-1} +O\left(v^{-\frac{4}{3}}\right)\right)\,,\quad 
\end{eqnarray}
where we denoted 
\[
v=\exp (3\theta-x)
\]
and $u_2$, $u_3$ are defined in terms of $p_i$ in (\ref{u_2_3}).

Since $U_i(x)$ constitute a complete set of solutions one can represent $\chi(x)$ as 
a linear combination
\begin{eqnarray}
\label{chiU}
\chi (x,\theta)=\sum_{n=1}^{3}Q_n(\theta)
\Gamma(p_{nj})\Gamma(p_{nk}) e^{-3p_n\theta}U_n(x,\theta),
\end{eqnarray}
where the important quantities $Q_n(\theta)$, based on general theory of linear 
differential equations, are expected to be entire functions of $\theta$ (and also of
parameters ${\bf p}$ dependence on which will be displayed explicitly only if necessary).  
The following, easy to check property plays an essential role in further discussion. 
Namely the Wronskian of any two solutions $f(x)$, $g(x)$ 
of the differential equation (\ref{diff_eq2}) 
\[
W[f(x),g(x)]\equiv f(x)g'(x)-g(x)f'(x)
\] 
satisfies the {\it adjoint} equation, i.e. the one obtained by reversing the signs $\bf{p}\rightarrow -\bf{p}$ 
and $\Lambda^3\rightarrow -\Lambda^3$. Taking inspiration from this property, it is then possible to show exactly that
\begin{eqnarray}
\label{Wr_chi_chi}
Wr\left[\chi (x,\theta+\frac{i \pi}{3}),\chi (x,\theta-\frac{i \pi}{3})\right]=
-\frac{i}{2 \pi }\bar{\chi}(x,\theta)\,,
\end{eqnarray}
where $\bar{\chi}(\theta)=\chi(\theta,-\mathbf{p})$. In fact, the property entails that the l.h.s. of (\ref{Wr_chi_chi}) satisfies the differential equation 
(\ref{diff_eq2}) with substitution $\bf{p}\rightarrow -\bf{p}$.  Besides, by using the identity\footnote{It can be proven, for instance, by expanding both sides in powers of $e^{-x}$.}
\begin{eqnarray}
&&Wr\left[e^{\frac{2-a-b}{3}\,x}\,_0F_2(a,b,-e^{-x})\,,e^{\frac{2a-b-1}{3}\,x}\,
_0F_2(2-a,1-a+b,-e^{-x})\right]\nonumber\\
&&=(a-1)e^{\frac{1+a-2b}{3}\,x}\,_0F_2(b,1-a+b,e^{-x})\,, 
\label{Wr_F_F}
\end{eqnarray}
it is not difficult to show the match of the $x\rightarrow-\infty$ asymptotics of both sides. Of course, the combination of these two statements implies the equality (\ref{Wr_chi_chi}) everywhere.

Let us investigate the $x\rightarrow \infty$ limit of (\ref{Wr_chi_chi}). Taking into account 
(\ref{chiU}) and using the identity (\ref{Wr_F_F}) (with $x$ substituted by $-x$), 
we obtain the functional relations
\begin{eqnarray}
\label{Qc_QQ}
\frac{\sin(\pi p_{jk})}{2i \pi^2}\, \bar{Q}_n(\theta)=
Q_j\left(\theta+\frac{i \pi}{3}\right)Q_k\left(\theta-\frac{i \pi}{3}\right)-
Q_j\left(\theta-\frac{i \pi}{3}\right)Q_k\left(\theta+\frac{i \pi}{3}\right),\qquad 
\end{eqnarray}
where again, the bar on $Q_n$ indicates the sign change $\bf{p}\rightarrow -\bf{p}$
\[
\bar{Q}_n(\theta,{\bf p})\equiv Q_n(\theta,-\bf{p})
\] 
and  $(n,j,k)$ is a permutations of $(1,2,3)$.

At the end of this section let us establish the $\theta\rightarrow-\infty $ asymptotics 
of $Q_k(\theta)$ and  $\bar{Q}_k(\theta)$. Obviously, in this case both (\ref{U_solutions})
and  (\ref{V_solutions}) are approximate solutions of (\ref{diff_eq2}) at $x\sim 0$. 
Thus, comparison of (\ref{chi}) with (\ref{chiU}) ensures that for $\theta\ll 0$
\begin{eqnarray} 
\label{Q_asymp}
Q_k(\theta)\sim \frac{\exp (-3\theta p_k)}{4 \pi^2}\,;\qquad 
\bar{Q}_k(\theta)\sim \frac{\exp (3\theta p_k)}{4 \pi^2}\,.
\end{eqnarray}
It is easy to see that above asymptotic behavior is fully  consistent
with functional relations (\ref{Qc_QQ}). 
\subsection{$SU(3)$ version of Baxter's $TQ$ relation}
The functional relations (\ref{Qc_QQ}) suggest the following $SU(3)$ analog 
of Baxter's $TQ$ equations:
\begin{eqnarray}
\label{T_Q}
&T(\theta)Q_j\left(\theta-\frac{\pi i}{6}\right)\bar{Q}_k\left(\theta+\frac{\pi i}{6}\right)  
=\hspace{9.6cm}\\
&Q_j(\theta-\frac{5\pi i}{6})\bar{Q}_k\left(\theta+\frac{\pi i}{6}\right)
+Q_j\left(\theta+\frac{\pi i}{2}\right)\bar{Q}_k\left(\theta-\frac{\pi i}{2}\right)+
Q_j\left(\theta-\frac{\pi i}{6}\right)\bar{Q}_k\left(\theta+\frac{5\pi i}{6}\right)\nonumber  
\end{eqnarray} 
for $j,k\in \{1,2,3\}$ with $j\neq k$. To uncover the essence of this construction, notice 
that for a fixed pair of indices $(i,j)$ (\ref{T_Q}) can be thought as definition 
of function $T(\theta)$ in terms of $Q$'s. Then the nontrivial question is ``do 
other choices of $(j,k)$ lead to the same $T$?" Fortunately, elementary algebraic manipulations 
with the help of (\ref{Qc_QQ}) ensure that the answer is positive. As mentioned 
earlier, $Q_i(\theta)$ are  entire functions. A thorough analysis shows that due to  
(\ref{Qc_QQ}) all potential poles of $T(\theta)$ have zero residue. Thus $T(\theta)$ 
is an entire function too.  Details on proofs of above two statements can be 
found in appendix \ref{TQ_proof}.
The Bethe ansatz equations can be represented as (see equality (\ref{T_entire}))
\begin{eqnarray} 
\frac{Q_j(\theta_\ell-\frac{2\pi i}{3})\bar{Q}_k\left(\theta_\ell+\frac{\pi i}{3}\right)}{
Q_j\left(\theta_\ell+\frac{2\pi i}{3}\right)\bar{Q}_k\left(\theta_\ell-\frac{\pi i}{3}\right)}=-1\,,
\end{eqnarray} 
where $\theta_\ell$ are the zeroes of $Q_j(\theta)$.  

Functional relations similar to (\ref{Qc_QQ}) and (\ref{T_Q}) emerge also in the context of ODE/IM for 'minimal' 2d CFT with extra spin $3$ current ($W_3$ symmetry) \cite{Dorey:1999pv}, \cite{Bazhanov:2001xm}. From there we can extrapolate that our case might correspond to the special choice 
of Virasoro central charge $c=98$ for Toda CFT. In fact, this value of the central charge lies outside the region discussed in above references. Nevertheless, it should be possible to derive the corresponding TBA equations: we leave this task for future publication. 
\section{Quantum periods and prepotential from Floquet monodromies and extension of Zamolodchikov's conjecture}
\label{a_cycles}
\subsection{The Floquet-Bloch monodromy matrix}
\label{sect_MM}
Consider the basis of solutions $f_1(x)$, $f_2(x)$, $f_3(x)$ of (\ref{diff_eq2}) with standard 
initial conditions 
($n,k\in\{1,2,3\}$)
\begin{eqnarray}
\label{BC}
\left.f^{(k-1)}_n(x)\right|_{x=0} =\delta_{k,n} \, .
\end{eqnarray}
Since the functions $f_n(x+2\pi i)$ are solutions too, we can
define the monodromy matrix $M_{k,n}$ as
\begin{eqnarray} 
f_n(x+2\pi i)=\sum_{k=1}^3 f_k(x)M_{k,n}
\end{eqnarray}
Clearly
\begin{eqnarray}
\label{M_matrix}
M_{k,n}=f_n^{(k-1)}(2\pi i)\nonumber \, .
\end{eqnarray}
The solutions (\ref{f_series}) with $a\in \{a_1,a_2,a_3\}$ have diagonal monodromies 
and can be represented as certain linear combinations of $f_n(x)$. 
In other words the eigenvalues of the monodromy matrix $M_{k,n}$ must be identified with $\exp (2\pi ia_k)$, with $k=1,2,3$:
\begin{eqnarray} 
\label{spec_M}
Spec(M_{k,n})=\{\exp (2\pi ia_1),\exp (2\pi ia_2),\exp (2\pi ia_3)\} \, .
\end{eqnarray} 

For any fixed values of parameters $\Lambda$, ${\bf p}$, it is easy to integrate 
numerically the differential equation (\ref{diff_eq2}) with boundary conditions 
(\ref{BC}), find the matrix $M_{k,n}$ and then its eigenvalues $\exp (2\pi ia_n)$. Taking into account Matone relation \cite{Matone:1995rx}, valid also in the presence of $\Omega$-background \cite{Flume:2004rp}, 
\begin{eqnarray}
u_2\equiv \langle \textbf{tr}\,\phi^2\rangle=\sum_{n=1}^3a_n^2+2q\partial_qF_{inst}(q,{\bf a}) \, ,
\end{eqnarray}   
we can access the deformed prepotential for any value of the coupling constant.
\subsection{Comparison of the instanton counting against numerical results}
Using formula of section \ref{NPF} it is straightforward to  calculate 
$\langle \textbf{tr}\,\phi^2\rangle$ or $\langle \textbf{tr}\,\phi^3\rangle$ as a 
power series in $q$.
Here are the 3-instanton results (it is assumed that $a_1+a_2+a_3=0 $ and 
by definition $a_{jk}\equiv a_j-a_k$) 
\begin{eqnarray} 
\label{tr_phi2}
\langle \textbf{tr}\,\phi^2\rangle &=&\sum_{k=1}^3a_k^2
-\frac{12(1-h_2)q}{\prod_{j<k}(a_{jk}^2-1)}
+\frac{P_{2,2}q^2}{\prod_{j<k}(a_{jk}^2-1)^3(a_{jk}^2-4)}+O(q)^4\qquad \\ 
\label{tr_phi3}
\langle \textbf{tr}\,\phi^3\rangle &=&\sum_{k=1}^3a_k^3
+\frac{54 h_3q}{\prod_{j<k}(a_{jk}^2-1)}
+\frac{P_{3,2}q^2}{\prod_{j<k}(a_{jk}^2-1)^3(a_{jk}^2-4)}\nonumber\\
&-&\frac{P_{3,3}q^3}{\prod_{j<k}(a_{jk}^2-1)^5(a_{jk}^2-4)(a_{jk}^2-9)}
+O(q)^4\,,
\end{eqnarray}  
where
\begin{eqnarray} 
h_2=\frac{a_1^2+a_2^2+a_3^2}{2}\,;\qquad h_3=-a_1a_2a_3\,,
\end{eqnarray}
and 
\begin{eqnarray} 
P_{2,2}=36 (220 - 1027 h_2 + 1659 h_2^2 - 698 h_2^3 - 958 h_2^4 + 1257 h_2^5 - 
521 h_2^6\hspace{3cm}\\
+ 68 h_2^7 - 13959 h_3^2 + 33804 h_2 h_3^2 - 
25434 h_2^2 h_3^2 + 5292 h_2^3 h_3^2\hspace{2.8cm}\nonumber\\
+ 297 h_2^4 h_3^2 + 13851 h_3^4 - 
5103 h_2 h_3^4)\hspace{3.5cm}\nonumber\\
P_{3,2}=-162 h_3 (455 - 2487 h_2 + 4602 h_2^2 - 3286 h_2^3 + 291 h_2^4 + 
573 h_2^5 - 148 h_2^6\hspace{2cm}\\
- 8073 h_3^2 + 14985 h_2 h_3^2 - 7695 h_2^2 h_3^2 + 
783 h_2^3 h_3^2 + 1458 h_3^4)\hspace{1.8cm}\nonumber\\
P_{3,3}=-108 h_3 (12078563 - 109310145 h_2 + 400164948 h_2^2 - 722480972 h_2^3
\hspace{2.5cm}\\ + 
538752402 h_2^4 + 275687658 h_2^5 - 946955868 h_2^6 + 
865056708 h_2^7\hspace{3.1cm}\nonumber\\ - 391259133 h_2^8 + 81882223 h_2^9 - 2063856 h_2^{10} - 
1715472 h_2^{11} + 162944 h_2^{12}\hspace{1.55cm} \nonumber\\- 984855213 h_3^2 + 
6130798389 h_2 h_3^2 - 14569978437 h_2^2 h_3^2+ 
16850898261 h_2^3 h_3^2\hspace{0.85cm} \nonumber\\ - 9439886367 h_2^4 h_3^2 + 
1593033399 h_2^5 h_3^2 + 730653777 h_2^6 h_3^2- 352792017 h_2^7 h_3^2 \hspace{1.1cm}\nonumber\\ + 
42690240 h_2^8 h_3^2 - 562032 h_2^9 h_3^2 + 7812512937 h_3^4 - 
22941081063 h_2 h_3^4 \hspace{2.1cm}\nonumber\\ + 24720233994 h_2^2 h_3^4 - 
11808597150 h_2^3 h_3^4 + 2295385533 h_2^4 h_3^4- 
64422459 h_2^5 h_3^4  \hspace{0.65cm}\nonumber\\- 14031792 h_2^6 h_3^4 - 3311723799 h_3^6 + 
3321565299 h_2 h_3^6 - 982634409 h_2^2 h_3^6 \hspace{1.7cm}\nonumber\\ + 65800269 h_2^3 h_3^6 + 
29760696 h_3^8) \hspace{2.6cm}\nonumber 
\end{eqnarray}
We have calculated also $4$ and $5$ instanton corrections, but the formulae are too lengthy 
to be presented here.
 
By means of numerical integration of the differential equation (\ref{diff_eq2}) along the 
line indicated in section \ref{sect_MM} we have computed the eigenvalues of monodromy matrix
(\ref{M_matrix}) for several values of the instanton parameter $q=\Lambda^6$, namely for 
the values 
\begin{eqnarray}
\label{Lambda_range}
\Lambda=\exp \left(\frac{k-1}{20}-5\right) \,,\qquad k=1,2,\cdots,120\,,
\end{eqnarray}  
and fixed values of parameters 
\[
p_1=0.12\,; \qquad p_2=0.17\,; \qquad p_3=-0.29\,.
\]
Due to identification (\ref{spec_M}) this allows to find the corresponding 
$A$-cycle periods $a_1,a_2,a_3$. In table \ref{tab:table1} we present some characteristic 
excerpt from the resulting data.  
\begin{table}[h!]
	\begin{center}
	
		\scalebox{0.9}{
\begin{tabular}{l|l|l}
\textbf{$\Lambda$} & $a_1$ & $a_2$\\
\hline
0.00822974704902 & 0.1200000000131 & 0.169999999982\\
0.0223707718562 & 0.1200000053049 & 0.169999992932 \\
0.0608100626252 & 0.1200021402877 & 0.169997148430\\
0.165298888222 & 0.1208841761521 & 0.168828966405\\
0.246596963942 & 0.1349151981823 & 0.151933010167 \\
0.272531793034 & 0.142136769453 - 0.019455438633 i & 0.142136769453 + 0.019455438633 i\\
0.449328964117 & 0.092117229441 - 0.135924390553 i & 0.092117229441 + 0.135924390553 i\\
0.740818220682 & 0.003727137475 - 0.568756791077 i & 0.003727137475 + 0.568756791077 i\\
1.22140275816 & 0.000899023180 - 1.071594057757 i & 0.000899023180 + 1.071594057757 i\\
2.01375270747 & 0.00036203460 - 1.78605985179 i & 0.00036203460 + 1.78605985179 i\\
3.32011692274 & 0.00013130957 - 2.96965962318 & 0.00013132399 + 2.96965962932 i
\end{tabular}}
\end{center}
\caption{The values $a_1$, $a_2$ obtained through numerical integration 
of the differential equation (\ref{diff_eq2}) with initial conditions (\ref{BC}) 
for $p_1=0.12$, $p_2=0.28$. }
\label{tab:table1}
\end{table}  

Inserting the values of $a_k$, $\Lambda$ in (\ref{tr_phi2}), (\ref{tr_phi3})
supplemented by $q^4$ and $q^5$ corrections we have calculated $\langle \textbf{tr}\,\phi^2\rangle$ and 
$\langle \textbf{tr}\,\phi^3\rangle$. The consistency requires that 
at small values of $q$ for at which instanton expansion is valid one should always obtain the same expectation values $\langle \textbf{tr}\,\phi^2\rangle 
=p_1^2+p_2^2+p_3^2=0.1274$ and $\langle \textbf{tr}\,\phi^3\rangle 
=p_1^3+p_2^3+p_3^3=-0.017748$. Table \ref{tab:table2} displays the results of actual
computations.   
\begin{table}[h!]
	\begin{center}
	
		\begin{tabular}{l|l|l}
			\textbf{$\Lambda$} & $\langle \textbf{tr}\,\phi^2\rangle$
			 & $\langle \textbf{tr}\,\phi^3\rangle$\\
			\hline
			0.00822974704902 & 0.1274000000000 & -0.0177480000000 \\
			0.0223707718562 & 0.1274000000000 & -0.0177480000000 \\
			0.0608100626252 & 0.1274000000000 & -0.0177480000000 \\
			0.165298888222 & 0.1274000000000 & -0.0177480000000 \\
			0.246596963942 & 0.1273999999998 & -0.0177480000000 \\
			0.272531793034 & 0.1273999999922 & -0.0177479999994 \\
			0.449328964117 & 0.1273774046391 & -0.0177462190257 \\
			0.740818220682 & 0.1313057536866 & -0.0178774876030
		\end{tabular}
	\end{center}
	\caption{The values $\langle \textbf{tr}\,\phi^2\rangle$, 
	$\langle \textbf{tr}\,\phi^3\rangle$  obtained by inserting  
	the values of $a_1$, $a_2$ from Table \ref{tab:table1} into (\ref{tr_phi2}), 
	(\ref{tr_phi3}) supplemented by $q^4$ and $q^5$ corrections. 
	To be compared with (by definition) $\langle \textbf{tr}\,\phi^2\rangle 
	=p_1^2+p_2^2+p_3^2=0.1274$ and $\langle \textbf{tr}\,\phi^3\rangle 
	=p_1^3+p_2^3+p_3^3=-0.017748$.}
	\label{tab:table2}
\end{table} 
One expects an essential deviation from the instanton series starting from 
the value of $\Lambda$ at which the polynomial 
\[
\left(z^3-\frac{u_2}{2}\,z-\frac{u_3}{3}\right)^2-4 \Lambda^6
\]
acquires coinciding zeros, i.e. at the point where its discriminant
\[
1024 \Lambda ^{18} \left(216 \Lambda ^6-72 \Lambda ^3 u_3-u_2^3+6 u_3^2\right) 
\left(216 \Lambda ^6+72 \Lambda ^3 u_3-u_2^3+6 u_3^2\right)
\]
vanishes. Such points correspond to massless dyons or monopoles. It is easy to 
check that within the range of $\Lambda$ 
(\ref{Lambda_range}) the only zero is at $\Lambda=0.1822359934629\cdots $ for which
the last factor of discriminant vanishes. And in fact, inspecting table \ref{tab:table2} 
one sees that for the greater values of $\Lambda$'s, the mismatch becomes significant 
while for smaller values the agreement is quite impressive. 
Notice also from table \ref{tab:table1} that for $\Lambda>0.24659696394 $ 
we encountered complex values for $a_{1}$ and $a_{2}$. 
\subsection{Extension of Zamolodcikov's conjecture to $SU(3)$}
The simpler case of the gauge group $SU(2)$ has been analyzed recently in 
\cite{Fioravanti:2019vxi}. In this case one has to deal with the Mathieu equation. Corresponding $TQ$ relation was investigated in \cite{Zamolodchikov:2000unpb}, where Al. Zamolodchikov conjectured (and demonstrated numerically) an elegant relationship 
between $T$-function and Floquet exponent $\nu$ of Mathieu equation:
\begin{eqnarray} 
T=\cosh (2\pi \nu) \, .
\end{eqnarray}
Here we suggest a natural extension of Zamolodchikov's conjecture for $SU(3)$ case:
\begin{eqnarray}
\label{conj_su3}
T(\theta)=\sum_{n=1}^3e^{2 \pi i a_n} \, .
\end{eqnarray} 
Notice, that at $\theta\ll 0$ the asymptotic (\ref{Q_asymp}) leads to
\begin{eqnarray} 
\label{T_asym}
T(\theta)\sim\sum_{n=1}^3e^{2 \pi i p_n}  \,\, ,
\end{eqnarray} 
which is consistent with (\ref{conj_su3}), since for $\theta\ll 0$ instanton 
corrections disappear and $a_k$ coincides with $p_k$.
\section{Few perspectives}
It would be very interesting to have a TBA for our case and check our 
conjecture (\ref{conj_su3}) as it was done by Al. Zamolodchikov in 
\cite{Zamolodchikov:2000unpb}. Actually, even relevant would be a gauge TBA that may shed light on the dual $B$-cycle periods 
$\bf{a_D}$ along the route presented in \cite{Fioravanti:2019vxi} for the $SU(2)$ case (see also the presentation of \cite{Gaiotto:2014bza} and \cite{Grassi:2019coc}).

It is well known that  pure ${\cal N}=2$ $SU(N)$ theories with $N>2$ are endowed with special points 
in there moduli spaces of vacua at which   mutually nonlocal dyons become massless (Argyres-Douglas points) \cite{ Argyres:1995jj}.
It would be interesting to investigate this phenomenon within our approach (for the NS regime of the $\Omega$ background).



Furthermore, for generic groups of gauge theories (starting with $SU(2)$ and general Liouville ODE/IM correspondence) it is very intriguing to investigate the form of 'potentials' of the ODE describing excited states of the IM ({\it cf.} \cite{BLZ-excited-ode-im, DF-excited-ode-im} for what we know about 'ordinary' ODE/IM). In fact, the latter should be obtainable also via analytic continuation in the parameters/moduli, and thus, these non-trivial monodromies would be of great interest in gauge theories.


Of course, it is very plausible that the imaginable generalizations of our results, and in particular of (\ref{conj_su3}), might hold for arbitrary $SU(N)$ gauge groups. In fact, for those higher order differential equations we shall have also the enlightening treatment of a mathematically similar problem with only one irregular singularity (at $\infty$), the case of gluon scattering amplitudes/Wilson loops at strong coupling in planar ${\cal N}= 4$ SYM \cite{TBA-amp1,TBA-amp3}. Because of its different physical nature, this problem allow for a beautiful and all-coupling exact Operator Product Expansion \cite{TBA-amp2,BSV1}, whose strong coupling limit reproduces interestingly the integrable TBA \cite{FPR} of \cite{TBA-amp1,TBA-amp3}: the similar mathematical structures and ideas of these two different fields should bear fruit, in future, for a deeper understanding.  


   


\section*{Acknowledgments}
HP and RP are grateful to G. Sarkissian and R. Mkrtchyan for many useful discussions.  
DF would like to thank D. Gregori, M. Rossi, R. Tateo for stimulating discussions. This work has been partially supported by the grants: research project 18T-1C340 from Armenian State Committee of Science, GAST (INFN), the MPNS-COST Action MP1210, the EC Network Gatis and the MIUR-PRIN contract 2017CC72MK$\textunderscore$003. D.F. thanks the GGI for Theoretical Physics for invitation to the workshop 'Supersymmetric Quantum Field Theories in the Non-perturbative Regime'.



\begin{appendix}
\section{Proving the $TQ$ relations}
\label{TQ_proof}
In this appendix we prove that different choices of indices $j\ne k $ in 
$TQ$ relations (\ref{T_Q}) are consistent with $QQ$ relations  (\ref{Qc_QQ}). 
 
For example let us choose $j=1$, $k=2$ 
\begin{eqnarray}
\label{T_Q_12}
&T(\theta)Q_1\left(\theta-\frac{\pi i}{6}\right)\bar{Q}_2\left(\theta+\frac{\pi i}{6}\right)  
=\\
&Q_1(\theta-\frac{5\pi i}{6})\bar{Q}_2\left(\theta+\frac{\pi i}{6}\right)
+Q_1\left(\theta+\frac{\pi i}{2}\right)\bar{Q}_2\left(\theta-\frac{\pi i}{2}\right)+
Q_1\left(\theta-\frac{\pi i}{6}\right)\bar{Q}_2\left(\theta+\frac{5\pi i}{6}\right),\nonumber  
\end{eqnarray} 
and  $j=1$, $k=3$  
\begin{eqnarray}
\label{T_Q_13}
&T(\theta)Q_1\left(\theta-\frac{\pi i}{6}\right)\bar{Q}_3\left(\theta+\frac{\pi i}{6}\right)  
=\\
&Q_1(\theta-\frac{5\pi i}{6})\bar{Q}_3\left(\theta+\frac{\pi i}{6}\right)
+Q_1\left(\theta+\frac{\pi i}{2}\right)\bar{Q}_3\left(\theta-\frac{\pi i}{2}\right)+
Q_1\left(\theta-\frac{\pi i}{6}\right)\bar{Q}_3\left(\theta+\frac{5\pi i}{6}\right),\nonumber  
\end{eqnarray} 
Multiplying (\ref{T_Q_12}) by $\bar{Q}_3\left(\theta+\frac{\pi i}{6}\right) $, 
(\ref{T_Q_13}) by $\bar{Q}_2\left(\theta+\frac{\pi i}{6}\right)$ 
and taking difference, the right hand side becomes
\begin{eqnarray}
&\bar{Q}_3\left(\theta+\frac{\pi i}{6}\right)\left(Q_1\left(\theta+\frac{\pi i}{2}\right)\bar{Q}_2\left(\theta-\frac{\pi i}{2}\right)+
Q_1\left(\theta-\frac{\pi i}{6}\right)\bar{Q}_2\left(\theta+\frac{5\pi i}{6}\right)\right)-
\nonumber\\
&\bar{Q}_2\left(\theta+\frac{\pi i}{6}\right)\left(
Q_1\left(\theta+\frac{\pi i}{2}\right)\bar{Q}_3\left(\theta-\frac{\pi i}{2}\right)+
Q_1\left(\theta-\frac{\pi i}{6}\right)\bar{Q}_3\left(\theta+\frac{5\pi i}{6}\right)
\right)=\\
&Q_1\left(\theta+\frac{\pi i}{2}\right)
\left(
\bar{Q}_3\left(\theta+\frac{\pi i}{6}\right)\bar{Q}_2\left(\theta-\frac{\pi i}{2}\right)-
\bar{Q}_2\left(\theta+\frac{\pi i}{6}\right)\bar{Q}_3\left(\theta-\frac{\pi i}{2}\right)\right)+
\nonumber\\
&Q_1\left(\theta-\frac{\pi i}{6}\right)\left(
\bar{Q}_3\left(\theta+\frac{\pi i}{6}\right)\bar{Q}_2\left(\theta+\frac{5\pi i}{6}\right)-
\bar{Q}_2\left(\theta+\frac{\pi i}{6}\right)\bar{Q}_3\left(\theta+\frac{5\pi i}{6}\right)
\right)\,.
\nonumber
\end{eqnarray}
Obviously the last expression vanishes due to (\ref{Qc_QQ}) as consistency requires.

Now let us show that $T(\theta)$ does not have any pole. 
It follows from (\ref{T_Q_12}), that a potential pole of $T(\theta)$ can be found 
either among the zeros of $Q_1\left(\theta-\frac{\pi i}{6}\right)$ or 
$\bar{Q}_2\left(\theta+\frac{\pi i}{6}\right)$. For definiteness let us assume that 
it belongs to zero set of  $Q_1\left(\theta-\frac{\pi i}{6}\right)$ (the other option 
can be considered in completely analogues manner). Let $Q_1(\theta_\ell)=0$. 
Then for $\theta=\theta_\ell+\frac{i\pi}{6}$ the r.h.s. of (\ref{T_Q_12}) is equal to 
\begin{eqnarray}
\label{T_entire}
&Q_1(\theta_\ell-\frac{2\pi i}{3})\bar{Q}_2\left(\theta_\ell+\frac{\pi i}{3}\right)
+Q_1\left(\theta_\ell+\frac{2\pi i}{3}\right)\bar{Q}_2\left(\theta_\ell-\frac{\pi i}{3}\right)=\\
&\frac{2i \pi^2}{\sin(\pi p_{31})}\left(Q_1(\theta_\ell-\frac{2\pi i}{3})Q_1\left(\theta_\ell+\frac{2\pi i}{3}\right)Q_3\left(\theta_\ell\right)-
Q_1\left(\theta_\ell+\frac{2\pi i}{3}\right)Q_1(\theta_\ell-\frac{2\pi i}{3})Q_3\left(\theta_\ell\right)\right)=0,
\nonumber
\end{eqnarray}
where the relations  (\ref{Qc_QQ}) for both 
$\bar{Q}_2\left(\theta_\ell\pm\frac{\pi i}{3}\right)$ has been used. So $T(\theta_\ell)$ is finite, 
thus proving that $T(\theta)$ is entire.
\end{appendix}
\bibliographystyle{JHEP}
\providecommand{\href}[2]{#2}
\begingroup\raggedright
