\section{I. Introduction}

  It is now recognized that cross-helicity
$$
H_{cr} = \int {\bf v}({\bf x},t)\cdot {\bf B}({\bf x},t) d{\bf r}  \eqno{(1)}
$$
plays important role in real MHD turbulence: in laboratory experiments, in atmosphere, in solar wind and in solar physics \cite{dmv}-\cite{pou1} (the cross-helicity is usually considered as a measure of relative contribution of the Alfv\'enic waves). Even for the cases with zero (or negligible) average cross-helicity
$$
h_{cr} = \langle  {\bf v}({\bf x},t)\cdot {\bf B}({\bf x},t) \rangle   \eqno{(2)}
$$
spatially localised cross-helicity density $({\bf v}({\bf x},t)\cdot {\bf B}({\bf x},t))$ can be rather large \cite{mene},\cite{matt} and, as it  will be shown below, can also play important role in the global MHD dynamics. Since the non-zero cross-helicity is naturally related to lack of the reflection symmetry (unlike the velocity field ${\bf v}$, which is a polar vector, the magnetic field ${\bf B}$ is an axial vector) the last phenomenon represents a kind of the spontaneous breaking of the reflection symmetry.\\

 The role of the cross-helicity is related to the fact that the cross-helicity is an invariant of the non-dissipative (ideal) MHD dynamics. It will be shown below that chaotic (coherent) nature of the MHD dynamics results in an extended inertial range dominated by the average cross-helicity Eq. (2) or by a second-order moment of cross-helicity fluctuations, which is also an invariant of the ideal MHD dynamics (MHD analogue of the Levich-Tsinober invariant \cite{lt},\cite{fl},\cite{l}) and can have a finite non-zero value even when the average cross-helicity is equal to zero. The extended inertial range has been described in the terms of distributed chaos approach and penetrates into near dissipation range of scales.

\section{II. Cross-helicity vs. magnetic energy}

   At the onset of turbulence in plasmas and  fluid dynamics deterministic chaos is often related the exponential power spectra \cite{mm},\cite{kds}
$$
E(k) \propto \exp-(k/k_c)  \eqno{(3)}
$$
where $k$ is wavenumber and $k_c=$ constant. Figure 1, for instance, shows (in the log-log scales) magnetic energy spectrum obtained in a recent direct numerical simulation (DNS) of the isotropic homogeneous MHD turbulence with a large-scale {\it deterministic} injection (forcing) of kinetic energy (see next section for more detail description). The spectral data used for the Fig. 1 were taken from Fig. 3 of the Ref. \cite{step}. The dashed curve indicates correspondence to the exponential spectrum Eq. (3) and the dotted arrow indicates position of the scale $k_c$. \\

\begin{figure} \vspace{-0.7cm}\centering
\epsfig{width=.45\textwidth,file=fig1.eps} \vspace{-3.8cm}
\caption{ Magnetic energy spectrum for the DNS with {\it deterministic} kinetic energy injection (forcing). }
\end{figure}

      For a random forcing, however, the parameter $k_c$ in the Eq. (3) fluctuates and in order to compute the spectra we need in ensemble averaging
$$
E(k) \propto \int_0^{\infty} P(k_c) \exp -(k/k_c)dk_c \propto \exp-(k/k_{\beta})^{\beta} \eqno{(4)}
$$
with certain probability distribution $P(k_c)$. The stretched exponential spectrum in the Eq. (4) is a natural generalization of the exponential spectrum Eq. (3). \\

  One can estimate the probability distribution $P(k_c)$ for large $k_c$ from the Eq. (4) \cite{jon}
$$
P(k_c) \propto k_c^{-1 + \beta/[2(1-\beta)]}~\exp(-\gamma k_c^{\beta/(1-\beta)}) \eqno{(5)}
$$

  Let us consider a scaling relationship between the parameter $k_c$ and characteristic magnetic field strength $B_c$
$$
B_c \propto  k_c^{\alpha}   \eqno{(6)}
$$

 In the case when $B_c$ has Gaussian distribution (with zero mean) a relationship between $\beta$ and $\alpha$
$$
\beta = \frac{2\alpha}{1+2\alpha}  \eqno{(7)}
$$
can be readily obtained from the Eqs. (5-6).

   In order to obtain value of $\alpha$ (and, consequently, value of $\beta$) let us use the dimensional considerations for the inertial range of scales. Namely
$$
B_c \propto    |\langle  {\bf v}\cdot {\bf B} \rangle|~\varepsilon^{-1/3}~ k_c^{1/3} \eqno{(8)}
$$
where $\varepsilon$ is the total energy (kinetic plus magnetic) dissipation rate (cf. the Corrsin-Obukhov approach for scaling of passive scalar spectrum Ref. \cite{my} and also Ref. \cite{bs}) . In the inertial range of scales the average cross-helicity $h_{cr} =  \langle  {\bf v}\cdot {\bf B} \rangle$ and the $\varepsilon$ can be considered as adiabatic invariants. Since the $\alpha = 1/3$ in this case we obtain from the Eq. (7) $\beta = 2/5$ i.e.
$$
E(k) \propto \exp-(k/k_{\beta})^{2/5}   \eqno{(9)}
$$

\section{III. Direct numerical simulations - I}

 Dynamics of an electrically conducting incompressible fluid can be described by the MHD equations
$$
 \frac{\partial {\bf v}}{\partial t} = - {\bf v} \cdot \nabla {\bf v}
    -\frac{1}{\rho} \nabla {\cal P} - [{\bf b} \times (\nabla \times {\bf b})] + \nu \nabla^2  {\bf v} + {\bf f_v} \eqno{(10)}
$$
$$
\frac{\partial {\bf b}}{\partial t} = \nabla \times ( {\bf v} \times
    {\bf b}) +\eta \nabla^2 {\bf b} +  {\bf f_b} \eqno{(11)}
$$
$$
\nabla \cdot {\bf v}=0, ~~~~~~~~~~~\nabla \cdot {\bf b}=0,  \eqno{(12)}
$$
where the velocity field ${\bf v}$ and normalized magnetic field ${\bf b} = {\bf B}/\sqrt{\mu_0\rho}$ have the same dimension (the so-called Alfv\'enic units), $ {\bf f}_v$ and ${\bf f}_b$ are forcing functions of the velocity and magnetic field respectively.
\begin{figure} \vspace{-2cm}\centering
\epsfig{width=.45\textwidth,file=fig2.eps} \vspace{-4.2cm}
\caption{ Magnetic energy spectrum for $C=0.3$. }
\end{figure}

\begin{figure} \vspace{-0.4cm}\centering
\epsfig{width=.45\textwidth,file=fig3.eps} \vspace{-4.3cm}
\caption{Magnetic energy spectrum for $C=0.6$.}
\end{figure}

  In the above mentioned DNS (Fig. 1) a statistically steady homogeneous and isotropic MHD-turbulence was simulated in a cubic volume using periodic boundary conditions. The {\it deterministic} large-scale forcing terms ${\bf f_b} =0$ and
$$
{\bf f_v} = \frac{ {\bf v}~}{|{\bf v}|^2},  \eqno{(13)}
$$
were applied in the wavenumber range $2 \leq k \leq 4$. The dissipative parameters were $\nu =0.008$ and $\eta = 0.004$.\\

   In the DNS reported in recent papers Refs. \cite{tit},\cite{mag} the deterministic forcing was replaced by a random one (with a controlled level of cross-helicity), which in the Fourier space can be written as
$$
{\bf f_v} ({\bf k}) = (\varepsilon_{s} (1-C)/2)^{1/2} {\bf e_v} ({\bf k}) + (\varepsilon_{s}C/2)^{1/2} {\bf e}_{cr} ({\bf k}),   \eqno{(14)}
$$
$$
{\bf f_b} ({\bf k}) = (\varepsilon_{s} (1-C)/2)^{1/2} {\bf e_b} ({\bf k}) + (\varepsilon_{s}C/2)^{1/2} {\bf e}_{cr} ({\bf k}),   \eqno{(15)}
$$
where
$$
{\bf e_{v,b,cr}}({\bf k}) = \frac{{\bf k}\times {\bf i_{v,b,cr}}}{|{\bf k}\times {\bf i_{v,b,cr}}|},   \eqno{(16)}
$$
${\bf i_{v,b,cr}}$ is random unit vector (unique for each subscript ${\bf v}, {\bf b}, {\bf cr}$) and updated whenever the forcing is initialized, $\varepsilon_{s} $ is total power of energy sources, $C=\varepsilon_{cr}/\varepsilon_{s}$ is level of injection of relative cross-helicity ($\varepsilon_{cr}$ is power of the sources of the cross-helicity). The large-scale forcing was applied in the wavenumber range $1 \leq k \leq 3$.

   Figures 2 and 3 show spectra of magnetic energy observed in the DNS reported in the Ref. \cite{tit} for $C=0.3$ and $C=0.6$ respectively (the spectral data were taken from the Fig. 2b of the Ref. \cite{tit}). The Reynolds and magnetic Reynolds numbers $Re =Re_m = 2094$. The dashed curves indicate correspondence to the stretched exponential spectrum Eq. (9).\\

   It is also interesting to consider a free (without forcing) decaying MHD turbulence with initially strong cross-correlation. Results of a DNS of such kind were reported in the Ref. \cite{baer}, where the so-called 3D Orszag-Tang flow \cite{spb} (with equal initial kinetic and magnetic energy and with the initial global cross-helicity normalized by total energy equal to 0.405) was taken as an initial condition.

     Figures 4 and 5 show spectra of magnetic energy observed in this DNS for the time of decay $t=1$ and $t=10$ (in the DNS terms) respectively. The spectral data were taken from Fig. 12 of the Ref. \cite{spb}. The dashed curves indicate correspondence to the stretched exponential spectrum Eq. (9).
\begin{figure} \vspace{-1.5cm}\centering
\epsfig{width=.45\textwidth,file=fig4.eps} \vspace{-4.3cm}
\caption{Magnetic energy spectrum for $t=1$.}
\end{figure}

\begin{figure} \vspace{-1.9cm}\centering
\epsfig{width=.45\textwidth,file=fig5.eps} \vspace{-3.9cm}
\caption{Magnetic energy spectrum for $t=10$.}
\end{figure}
\begin{figure} \vspace{-0.4cm}\centering
\epsfig{width=.45\textwidth,file=fig6.eps} \vspace{-3.9cm}
\caption{Magnetic energy spectrum in the Earth's magnetosheath.}
\end{figure}

\section{IV. Cross-helicity in the Earth's Magnetosheath}

  The magnetohydrodynamic description is considered as an adequate one for the large-scale processes in the solar wind and in the Earth's magnetospheric region located downstream of the bow shock (the so-called Earth's magnetosheath). In the recent Ref. \cite{ban} the magnetic energy spectrum was computed using data obtained by the Magnetospheric Multiscale  (MMS) Mission spacecraft operated in the Earth's magnetosheath \cite{bur}. The low Mach number and density fluctuations (Table 1 of the Ref. \cite{ban}) justify applicability of the incompressible magnetohydrodynamics for the large scales. The estimated by the authors normalized cross-helicity $\sigma_c \simeq 0.24$ (Table 2 of the Ref. \cite{ban}) was smaller than those which can be observed in the solar wind but was certainly not negligible (let us recall that maximal value of $\sigma_c =1$).

    Figure 6 shows the magnetic energy spectrum against frequency (the spectral data were taken from the Fig. 2 of the Ref. \cite{ban}). Since in this case the wind's mean velocity $|\langle {\bf V} \rangle|$ in the spacecraft frame is considerably larger than the representative velocity fluctuations the Taylor "frozen" hypothesis $k \simeq 2\pi f/ |\langle {\bf V}| \rangle$ can be applied (see, for instance, Refs. \cite{hb},\cite{tbn} and references therein). Following to this hypothesis the temporal dynamics measured by the probe merely reflects the spatial one convected past the probe by the mean (nearly constant) velocity. Therefore, it is not a true frequency spectrum but actually a wavenumber spectrum (about intrinsic time dependence of the magnetic field see the Section VIII). Hence, the dashed curve in the Fig. 6 indicates correspondence to the stretched exponential {\it wavenumber} spectrum Eq. (9) and the scale $k_{\beta}\simeq 2\pi f_{\beta}/ |\langle {\bf V} \rangle|$ (it is clear that the $k_{\beta}$ corresponds to large-scale spatial structures in this case). The vertical black bar indicates the ion inertial length $d_i$ ($kd_i=1$).

\section{V. Cross-helicity vs. total energy}

   It was already mentioned in the Introduction that since the non-zero cross-helicity is related to lack of the reflection symmetry (unlike the velocity field ${\bf v}$, which is a polar vector, the magnetic field ${\bf B}$ is an axial vector) the $H_{cr}$ and $h_{cr}$ Eqs. (1-2) can be non-zero only when the global reflection symmetry is broken. Their role in MHD is similar to the role of the ideal invariants based on the hydrodynamic helicity $({\bf v} \cdot  {\boldsymbol \omega})$ (where ${\boldsymbol \omega} = \nabla \times {\bf v}$ is vorticity) in ordinary hydrodynamics (see the Refs. \cite{lt},\cite{b1} and for reviews the Refs. \cite{l},\cite{mt} and references therein). In particular, in the case of global reflection symmetry $H_{cr}=h_{cr}=0$. However, the higher (even) moments of the cross-helicity fluctuations can nevertheless be finite and constant, due to spatially localized lack of the reflection symmetry (with mutual compensation of contribution of the spatial areas with different sign of the cross-helicity distribution into the global cross-helicty). Therefore, the higher moments are of a special interest. Moreover, even in
the cases of lack of the global reflection symmetry (when $H_{cr}$ and $h_{cr}$ are non-zero) the second moment can dominate extended inertial range of the total energy spectrum (as it will be shown below, cf. also the Ref. \cite{b1}).  \\
\begin{figure} \vspace{-1.7cm}\centering
\epsfig{width=.45\textwidth,file=fig7.eps} \vspace{-4.1cm}
\caption{Total energy spectrum for $C=0.3$.}
\end{figure}
\begin{figure} \vspace{-1.85cm}\centering
\epsfig{width=.44\textwidth,file=fig8.eps} \vspace{-3.7cm}
\caption{Total energy spectrum for $C=0.6$.}
\end{figure}

  Since dynamics of the magnetic field ${\bf b}$ in the non-dissipative case
$$
\frac{\partial {\bf b}}{\partial t} = \nabla \times ( {\bf v} \times
    {\bf b})   \eqno{(17)}
$$
has the same form as for vorticity ${\boldsymbol \omega}$ (the initial conditions for this equation, of course, belong to a much wider class than those allowed for the vorticity's dynamical equation \cite{mt}) the same consideration that was applied to the helicity moments \cite{lt},\cite{mt} can be applied to the cross-helicity moments as well. Namely, let us divide the entire volume of motion into cells $V_j$ with boundary conditions ${\bf b} \cdot {\bf n}=0$ on the bounding surfaces $S_j$ moving with the fluid. Then the 'localized' in the volume $V_j$ cross-helicity
$$
H_{cr,j} = \int_{V_j} {\bf v}({\bf x},t)\cdot {\bf b}({\bf x},t) d{\bf r}  \eqno{(18)}
$$
is a non-dissipative (ideal) invariant of the motion. The second order moment of the cross-helicity distribution in this case can be defined as
$$
I_{cr} = \lim_{V \rightarrow  \infty} \frac{1}{V} \sum_j H_{cr,j}^2  \eqno{(19)}
$$
and it is, due to its construction, a non-dissipative (ideal) invariant of the motion (it is interesting to compare this approach to the multifractal one, see for instance Ref. \cite{bt} and references therein). \\

  Let us consider a characteristic 'velocity' - $\tilde{v}_c$  defined for the total energy:
$$
{\mathcal E} = \frac{1}{2}({\bf v}^2 + {\bf b}^2)     \eqno{(20})
$$
(in the Alfv\'enic units the normalized magnetic field ${\bf b} = {\bf B}/\sqrt{\mu_0\rho}$ have the same dimension as velocity). From the dimensional considerations we obtain
$$
\tilde{v}_c \propto I_{cr}^{1/4} k_c^{3/4},  \eqno{(21)}
$$
i.e. $\alpha =3/4$ in this case (cf. Eq. (6)). Then from the Eq. (7) we obtain $\beta =3/5$, i.e. the total energy spectrum
$$
E(k) \propto \exp-(k/k_{\beta})^{3/5}.   \eqno{(22)}
$$

\section{VI. Direct numerical simulations - II}

   Let us again start from the direct numerical simulation described in the Ref. \cite{tit} (see the Section: "Direct numerical simulations - I" above, Figs. (2-3)), but now we will consider results obtained for the total energy spectrum and reported in recent Ref. \cite{mag}. The {\it random} large-scale forcing terms Eqs. (14-15) were applied in the wavenumber range $1 \leq k \leq 3$ and the Kolmogorov (dissipation) wavenumber $k_K = k_m \simeq 100$ ($Re = Re_m \simeq 2094$). Figures 7 and 8 show spectra of total energy for $C=0.3$ and $C=0.6$ respectively (the spectra correspond to the Fig. 1a of the Ref. \cite{mag}). The dashed curves indicate correspondence to the stretched exponential spectrum Eq. (22) and the dotted arrows indicate position of the scale $k_{\beta}$.

   One can see that the extended inertial range penetrates rather deep into the near dissipation range. This phenomenon can be related to the coherency with the large scales introduced by domination of the distributed chaos (cf. position of the scale $k_{\beta}$ in the Figs. (7-8)). It should be noted in this respect that the cross helicity is known as a factor that introduce a nonlocal interaction between the small- and large-scale fluctuations. The coherency results in the adiabatic invariance of the $I_{cr}$ up to the the dissipation scale (cf. also Refs. \cite{lst},\cite{mof}).  \\

   Let us also consider dynamics of the total energy in the freely decaying MHD turbulence with initially strong cross-correlation generated by the initial conditions taken in the form of the 3D Orszag-Tang flow (cf. the Section: "Direct numerical simulations - I" above, Figs. (4-5)). Figures 9 and 10 show spectra of total energy for the Taylor-Reynolds numbers $Re_{\lambda} = 45$ and  $Re_{\lambda} = 30$ respectively (the spectral data were taken from the Fig. 2a of the Ref. \cite{gibbon} and correspond to the early times of the decay, the $Re_{\lambda}$ is decreasing with time of decay). The magnetic Prandtl number $P_m$ was again taken equal to 1, i.e $\nu=\eta$. The dashed curves indicate correspondence to the stretched exponential spectrum Eq. (22) and the dotted arrows indicate position of the scale $k_{\beta}$. \\

   Finally, let us consider results of DNS with external/mean magnetic field reported in Ref. \cite{cho} (this case is  important for astrophysics). The imposed magnetic field was a moderate one $b_0 =1$ (in the terms of the DNS \cite{cho}). However, the statistically stationary MHD turbulence is anisotropic. The MHD turbulence was incompressible and was described by the Eq. (10-12) with $\nu=\eta$ and random isotropic forcing ${\bf f_v}$ (whereas ${\bf f_b} =0$). The large-scale random energy injection has its peak at $k \simeq 2.5$.
\begin{figure} \vspace{-1.5cm}\centering
\epsfig{width=.42\textwidth,file=fig9.eps} \vspace{-4cm}
\caption{Total energy spectrum for $Re_{\lambda} = 45$.}
\end{figure}
\begin{figure} \vspace{-0.5cm}\centering
\epsfig{width=.45\textwidth,file=fig10.eps} \vspace{-4.2cm}
\caption{Total energy spectrum for $Re_{\lambda} = 30$.}
\end{figure}
\begin{figure} \vspace{-1.5cm}\centering
\epsfig{width=.45\textwidth,file=fig11.eps} \vspace{-4.2cm}
\caption{Total energy spectrum for $b_0 =1$.}
\end{figure}
\begin{figure} \vspace{-0.5cm}\centering
\epsfig{width=.445\textwidth,file=fig12.eps} \vspace{-3.85cm}
\caption{Total energy spectrum for $C = 0$.}
\end{figure}
\begin{figure} \vspace{-0.5cm}\centering
\epsfig{width=.445\textwidth,file=fig13.eps} \vspace{-4.5cm}
\caption{Magnetic energy spectrum for $C = 0$.}
\end{figure}
   Figure 11 shows spectrum of total energy observed in this DNS for $b_0 =1$ (the spectral data were taken from the Fig. 16 of the Ref. \cite{cho}). The dashed curve indicates correspondence to the stretched exponential spectrum Eq. (22) and the dotted arrows indicate position of the scale $k_{\beta}$.

\section{VII. Spontaneous breaking of local reflection symmetry}

The cross-helicity density changes its sign at mirror reflection of the coordinate system (a pseudoscalar). Therefore, the global quantities  Eqs. (1-2) are non-zero only if there is a global lack of reflection (mirror) symmetry, and the non-zero average cross-helicity indicates the breakage of global reflection (mirror) symmetry.  However, even when the average cross-helicity is negligible the magnitude of its density can be locally large (in the vicinity of current and vorticity sheets, for instance, see the Rerfs. \cite{mene},\cite{matt2} and references therein). After the global averaging over the localized patches with negative and positive cross-helicity one has the average cross-helicity close to zero (there exits a view that the MHD turbulence with close to zero average cross-helicity can be generally considered as a superposition of these patches, see for instance Ref. \cite{pb} and references therein). Moreover, the patches can possess a kind of hierarchical structure. Namely, inside the patches there can exist smaller ones with more strong chirality (of different signs) and so on \cite{pb}. The above considered (adiabatic) invariant $I_{cr}$ Eq. (19) is rather adequate tool for description of this situation. In particular, as it was already mentioned, this (adiabatic) invariant can be finite even in the cases when the average cross-helicty is close to zero. \\

  To check applicability of this consideration let us return to the results obtained in the direct numerical simulation reported in the Ref. \cite{mag} (see previous Section) but now for $C = 0$, i.e. for the case with negligible average cross-helicity. Figure 12 shows spectrum of total energy obtained for $C=0$ (the spectrum corresponds to the Fig. 1a of the Ref. \cite{mag}). The dashed curve indicates correspondence to the stretched exponential spectrum Eq. (22) and the dotted arrow indicates position of the scale $k_{\beta}$. \\

  It is instructive to compare Fig. 12 with the Figs. (7-8) showing the total energy spectra for $C=0.3$ and $C=0.6$ (corresponding to considerable average cross-helicity).  This comparison allows us to conclude that the effect of the spontaneous breaking of the local reflection (mirror) symmetry indeed takes place in this case and can be described in the terms of the (adiabatic) invariant $I_{cr}$ (see also below).

    It should be noted that in the Alfv\'enic units the normalized magnetic field ${\bf b} = {\bf B}/\sqrt{\mu_0\rho}$ has the same dimension as velocity and the Eq. (21) can be replaced by equation
$$
b_c \propto I_{cr}^{1/4} k_c^{3/4}  \eqno{(23)}
$$
Therefore, unlike the cases with considerable average cross-helicty (see Figs. 2 and 3), for $C=0$ the {\it magnetic} energy spectrum has behaviour similar to that shown in the Fig. 12. Figure 13 shows the magnetic energy spectrum obtained for $C=0$ (the spectral data for the Fig. 13 were taken from the Fig. 2b of the Ref. \cite{tit}, as for the Figs. 2 and 3). The dashed curve indicates correspondence to the stretched exponential spectrum Eq. (22), whereas for the Figs. 2 and 3 the magnetic energy spectrum corresponds to the Eq. (9). \\

    Generally the MHD turbulence exhibits a wide variety of energy spectra depending on initial-boundary conditions and type of forcing. At present time the necessary and sufficient conditions for the spontaneous breaking of local reflection symmetry in MHD turbulence are not known.

\section{VIII. Spatio-temporal distributed chaos}

  Is the above considered distributed MHD chaos exclusively spatial or is it a spatio-temporal one? To answer this question let us replace the spatial (wavenumber)  Eq. (21) by its temporal (frequency) analogue using the dimensional considerations:
$$
\tilde{v}_c \propto I_{cr}^{1/7} f_c^{3/7}  \eqno{(24)}
$$
where $f_c$ is a characteristic frequency. It follows from the Eq. (24) that $\alpha =3/7$. Then from the Eq. (7) we obtain $\beta =6/13$, i.e. the frequency total energy spectrum is
$$
E(f) \propto \exp-(f/f_{\beta})^{6/13}.   \eqno{(25)}
$$
in this case.

   In the Alfv\'enic units the normalized magnetic field ${\bf b} = {\bf B}/\sqrt{\mu_0\rho}$ has the same dimension as velocity and the Eq. (24) can be replaced by equation
$$
b_c \propto I_{cr}^{1/7} f_c^{3/7},  \eqno{(26)}
$$
i.e. in this case the magnetic energy spectrum has the same form as the total energy spectrum Eq. (25). \\

  In Ref. \cite{dm} results of a direct numerical simulation of the isotropic homogeneous MHD turbulence of incompressible conducting fluid were reported and an Eulerian {\it frequency} spectrum of magnetic energy was constructed using 64 point-like probes set in a middle plane of 3D spatial box in a regular  $8\times8$ array (the Reynolds and magnetic Reynolds numbers were $Re=Re_m=400$). The long time series obtained with these probes were used to compute the average (over all probes) frequency spectrum of magnetic energy fluctuations. The initial state for this simulation was a random phased set of velocity and magnetic filed fluctuations (in equipartition) in the shell $1 \leq |{\bf k}| \leq 4$. The initial average cross-helicity was negligible in comparison with the energy. Due to the {\it uncorrelated} character of the random forcing ${\bf f_v}$ and ${\bf f_b}$ (in the range $1 \leq  k \leq 2$) there was no statistical injection of cross-helicity and the average cross-helicity was negligible also in the steady state of the MHD turbulence. Therefore, one can expect that the spontaneous breaking of local reflection symmetry can take place at this simulation. \\

\begin{figure} \vspace{-1.7cm}\centering
\epsfig{width=.445\textwidth,file=fig14.eps} \vspace{-4.1cm}
\caption{Frequency spectrum of magnetic energy.}
\end{figure}
   Figure 14 shows frequency spectrum of magnetic energy obtained in this DNS for the incompressible MHD turbulence without background magnetic field. The spectrum corresponds to the Fig. 2 (upper panel) of the Ref. \cite{dm}. The dashed curve indicates correspondence to the stretched exponential spectrum Eq. (25) and the dotted arrow indicates position of the scale $f_{\beta}$. One can see that the four decades range of the frequency spectrum supports the spatio-temporal character of the distributed MHD chaos and indicates the spontaneous breaking of the local reflection symmetry in this case (cf. Fig. 13 corresponding to the wavenumber spectrum).

\section{IX. Effects of mean magnetic filed}

  We have already mentioned results of DNS with a moderate imposed magnetic field (see Fig. 11). Presence of a strong mean magnetic field can result in fundamental changes in the properties of MHD turbulence. Even when one ignore anisotropy introduced by the mean magnetic field (suggesting, for instance, the local isotropy), as it was made in the pioneering papers \cite{ir},\cite{kr}, the dimensional parameter governing the inertial range of scales in the Kolmogorov theory - $\varepsilon$, should be replaced by the parameter $\varepsilon b_0$ (where $\varepsilon$ is the energy dissipation rate and the normalised mean magnetic field $b_0 = B_0/\sqrt{\mu_0\rho}$ has the same dimension as velocity). A vigorous discussion about validity of this replacement in solar wind, for instance, still takes place in interpretation of scaling properties of the modern spacecraft data (spectra $E(k) \propto k^{-5/3}$ versus spectra $E(k) \propto k^{-3/2}$) and a few new theoretical approaches were suggested (let us mention the Refs. \cite{gs},\cite{bold}).\\

\begin{figure} \vspace{-1.7cm}\centering
\epsfig{width=.445\textwidth,file=fig15.eps} \vspace{-4.1cm}
\caption{Magnetic energy spectrum in the MHD wind (plasma) tunnel.}
\end{figure}
  With the replacement $\varepsilon \rightarrow \varepsilon b_0$ the Eq. (8) should be replaced by the equation
$$
B_c \propto    |\langle  {\bf v}\cdot {\bf B} \rangle|~(\varepsilon b_0)^{-1/4}~ k_c^{1/4} \eqno{(27)}
$$
that gives $\alpha =1/4$ and, correspondingly, one obtains from the  Eq. (7) $\beta = 1/3$. Hence the magnetic energy spectrum is
$$
E(k) \propto \exp-(k/k_{\beta})^{1/3}   \eqno{(28)}
$$
in this case. \\

   Recent paper Ref. \cite{sbl} reports results obtained in a laboratory experiment with magnetic turbulent plasma in a MHD wind tunnel. The experiment was especially designed in order to model solar wind magnetohydrodynamics. Figure 15 shows magnetic energy spectrum (ensemble averaged) measured in this experiment (the spectral data for this figure were taken from Fig. 15a of the Ref. \cite{sbl}). The Taylor hypothesis (see Section IV) relates the frequency spectrum to analogous spatial (wavenumber) spectrum $E(k)$ and the dashed curve indicates correspondence to the stretched exponential spectrum Eq. (28).\\

      As for the solar wind itself the spacecraft Ulysses and Helios-1 measurements, made at high and low heliolatitudes respectively, provide a general picture of its magnetohydrodynamics from the distances 4.5 AU (Ulysses) to 0.3 AU  (Helios-1) from the Sun. These measurements in the high-speed streams show a strong similarity for the Ulysses and Helios-1 data (see, for instance, Ref. \cite{hb} and references therein).\\

      Figure 16 shows an example of magnetic energy spectra obtained by magnetometers of the Helios-1. The spectral data for the Fig. 16 were taken from Fig. 3 of the Ref. \cite{tbn} (full speed mapping). The dashed curve indicates correspondence to the stretched exponential spectrum Eq. (28). The dotted arrow indicates position of the scale $k_{\beta}$. \\

  Figures 17 and 18 show magnetic energy spectrum corresponding to the Ulysses data for the period 1993-1996yy at high solar wind speed and at high heliolatitudes. The spectral data for these figures were taken from Fig. 3 of the Ref. \cite{bran} for $1.5 < R < 2.8$ AU and for $2.8 < R < 4.5$ AU, respectively ($R$ is distance between the spacecraft and the Sun). Before averaging over the data sets the spectra were rescaled by factor $4\pi R^2$. The dashed curves correspond to the Eq. (28).

\begin{figure} \vspace{-3.3cm}\centering
\epsfig{width=.445\textwidth,file=fig16.eps} \vspace{-2.5cm}
\caption{Magnetic energy spectrum in the solar wind (the Helios-1 data).}
\end{figure}

\begin{figure} \vspace{-0.3cm}\centering
\epsfig{width=.445\textwidth,file=fig17.eps} \vspace{-4.1cm}
\caption{Magnetic energy spectrum in the solar wind (the Ulysses data for $1.5 < R < 2.8$ AU).}
\end{figure}

\begin{figure} \vspace{-1.5cm}\centering
\epsfig{width=.445\textwidth,file=fig18.eps} \vspace{-4.3cm}
\caption{Magnetic energy spectrum in the solar wind (the Ulysses data for $2.8 < R < 4.5$ AU).}
\end{figure}

\section{Acknowledgement}

I thank  E. Levich for stimulating discussions, A. Beresnyak, V. Pipin, D.D. Sokoloff and J. Shebalin for comments,
R. Grappin, R. Stepanov and V. Titov for sharing their data and additional information.


\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