\section{#1}\setcounter{equation}{0}}
\renewcommand{\theequation}{\thesection.\arabic{equation}}
\newcommand\ef[1]{(\ref{#1})}

\newdimen\tableauside\tableauside=1.0ex
\newdimen\tableaurule\tableaurule=0.4pt
\newdimen\tableaustep
\def\phantomhrule#1{\hbox{\vbox to0pt{\hrule height\tableaurule width#1\vss}}}
\def\phantomvrule#1{\vbox{\hbox to0pt{\vrule width\tableaurule height#1\hss}}}
\def\sqr{\vbox{%
  \phantomhrule\tableaustep
  \hbox{\phantomvrule\tableaustep\kern\tableaustep\phantomvrule\tableaustep}%
  \hbox{\vbox{\phantomhrule\tableauside}\kern-\tableaurule}}}
\def\squares#1{\hbox{\count0=#1\noindent\loop\sqr
  \advance\count0 by-1 \ifnum\count0>0\repeat}}
\def\tableau#1{\vcenter{\offinterlineskip
  \tableaustep=\tableauside\advance\tableaustep by-\tableaurule
  \kern\normallineskip\hbox
    {\kern\normallineskip\vbox
      {\gettableau#1 0 }%
     \kern\normallineskip\kern\tableaurule}%
  \kern\normallineskip\kern\tableaurule}}
\def\gettableau#1{\ifnum#1=0\let\next=\null\else
\squares{#1}\let\next=\gettableau\fi\next}

\tableauside=1.0ex
\tableaurule=0.4pt
\def\mathbb{E}{\mathbb{E}}

\def\mmm#1{$\spadesuit\spadesuit$ \textcolor{red}{MM: \emph{#1}}
  $\spadesuit\spadesuit$}

\newcommand{\figref}[1]{Fig.~\protect\ref{#1}}
\usepackage{mathtools}

\title{\Huge{\boldmath Instantons, renormalons and the theta angle in integrable sigma models}}

\author{Marcos Mari\~no,}
\author{Ramon Miravitllas}
\author{and Tom\'as Reis}

\affiliation{D\'epartement de Physique Th\'eorique et Section de Math\'ematiques\\
Universit\'e de Gen\`eve, Gen\`eve, CH-1211 Switzerland}

\emailAdd{Marcos.Marino@unige.ch}
\emailAdd{Ramon.MiravitllasMas@unige.ch}
\emailAdd{Tomas.Reis@unige.ch}

\abstract{Some sigma models which admit a theta angle are integrable at both $\vartheta=0$ and $\vartheta=\pi$. This includes the
well-known $O(3)$ sigma model and two families of coset sigma models studied by Fendley. We consider the ground state energy of
these models in the presence of a magnetic field, which can be computed with the Bethe ansatz. We obtain explicit results for its
non-perturbative corrections and we study the effect of the theta angle on them. We show that imaginary, exponentially small
corrections due to renormalons remain unchanged, while instanton corrections change sign, as expected. We find in addition corrections due to renormalons which also change sign as we turn on the theta angle.
Based on these results we present an explicit non-perturbative
formula for the topological susceptibility of the $O(3)$ sigma model
in the presence of a magnetic field, in the weak coupling limit.}

\begin{document}
\maketitle
\flushbottom

\sectiono{Introduction}

Perturbative series in quantum field theory have typically non-perturbative
corrections which are exponentially small in the coupling
constant. We expect these corrections to be due either to instantons or to renormalons. Instanton corrections
come from non-trivial saddle-points of the path integral and they are in principle accessible with
semiclassical methods. In practice, however, instanton calculus is plagued with all kinds of difficulties, the most important one being
perhaps the infrared (IR) divergences appearing in asymptotically free theories. Renormalons are even more
puzzling, since they have no known description in terms of saddle points of the path integral, and their current understanding
relies heavily on estimates from perturbation theory at large order.

One possible avenue for a better understanding of non-perturbative corrections is to
look at solvable models where one can extract them from the exact solution. The results
obtained in this way give a signpost to test non-perturbative methods.
One example of this strategy is the Seiberg--Witten solution of ${\cal N}=2$ supersymmetric Yang--Mills theory \cite{sw},
which provides a prediction for multi-instanton contributions in this theory. This prediction
was first addressed by conventional methods (see \cite{hollo-insts} for a review), and the effort to understand instantons
in this theory from first principles culminated in Nekrasov's work \cite{n}.
However, this is a supersymmetric success story in which many generic aspects
of non-perturbative effects can not be detected. In particular, ${\cal N}=2$ super Yang--Mills theory has
no renormalons, and instanton corrections are one-loop exact.

It has been known for a long time that certain two-dimensional, asymptotically free theories are integrable,
and their exact $S$-matrices have been conjectured. These include the $O(N)$ supersymmetric sigma model,
the Gross--Neveu model, and the principal chiral field (PCF). These models are expected to have renormalons,
and some of them admit instantons and other classical solutions. Therefore, they provide a rich laboratory where
one might obtain explicit results for non-perturbative effects, like instanton and renormalon corrections.
This issue has been addressed in previous works, by working at large $N$ \cite{fkw1,fkw2,dpmss} or by using numerical techniques
\cite{volin, volin-thesis,mr-ren,abbh1,abbh2}\footnote{There have been important works on non-perturbative corrections in two-dimensional quantum
field theories which do not rely on integrability, but only on large $N$ solvability. These include David's
pioneering work on renormalons \cite{david2,david3}, as well as \cite{itep2d,beneke-braun,shifman}. Another, more recent line of work
considers twisted compactifications of two-dimensional models, see e.g. \cite{dunne-unsal-cpn, cddu}.
However, the calculation of perturbative series and non-perturbative effects in \cite{dunne-unsal-cpn, cddu} and follow-up papers is based on an
uncontrolled truncation to an effective quantum mechanical model, and provides at best approximate results.}. In our previous
paper \cite{mmr-an} we have developed a general method to obtain {\it analytic, exact} results at finite $N$ for
exponentially small corrections to the free energy in these integrable
models\footnote{The recent work \cite{bbv} obtains analytic results for non-perturbative effects in the $O(4)$ sigma model.}.
Some of these corrections were identified with renormalon effects, and we showed in particular that their structure is more general
than the one put forward in the pioneering studies of Parisi \cite{parisi1,parisi2} and 't Hooft \cite{thooft}.
However, it was pointed out in \cite{mmr-an} that some of the corrections found in the $O(N)$
sigma model and in the PCF seem to correspond to instantons.
In particular, as emphasized in \cite{bbh}, the leading, real exponentially small correction in the $O(3)$ non-linear sigma model might
be due to the well-known instanton configurations in that theory \cite{bel-pol}.

One important and genuine non-perturbative effect in quantum field theory is the dependence on the
topological theta angle, which is completely invisible in perturbation theory. In some cases this dependence
is due to instantons, although large $N$ considerations \cite{witten-largeN} and explicit lattice
calculations \cite{giusti} in pure Yang--Mills theory show that the
physics of the theta angle can not be always explained by semiclassical methods.
Some integrable two-dimensional models, like
the $O(3)$ sigma model, admit also a theta angle. Although turning on a
generic theta angle is believed to
break integrability, it was observed in \cite{fz,zz-theta}
that the $O(3)$ sigma model with $\vartheta =\pi$ is still integrable. This picture was generalized by Fendley \cite{fendley} to
two families of coset sigma models: the $SU(N)/SO(N)$ models, and the $O(2P)/O(P) \times O(P)$ models.
In these theories there is a ${\mathbb Z}_2$ instanton number, and the theta angle can only take the values $\vartheta=0, \pi$.
Fendley conjectured that these models are also integrable for both values of $\vartheta$, and proposed exact expressions for their $S$-matrices.

In view of the integrability of all these models, one could try to use their exact solutions at $\vartheta=0,\pi$ to study the behavior of
non-perturbative effects as we change the theta angle. This issue was already addressed
in \cite{foz,fendley} by considering deformed versions of these models
(in the case of the $O(3)$ sigma model, this deformation is known as the ``sausage" model).
The calculation of non-perturbative effects in these deformed models is much simpler than in the
original ones, and this was part of the motivation to introduce the deformation in the first place.
However, the deformed models have various drawbacks: their spacetime physics is not transparent,
and they do not seem to have renormalon effects. In addition,
the original models are recovered from the deformed ones in a singular limit. This means in particular that the
free energy of the original, undeformed models cannot be obtained from the free energy of the deformed
models calculated in \cite{foz,fendley}. For this reason, the results of \cite{foz, fendley} do not give direct access to
non-perturbative corrections in the conventional sigma models.

In this paper we study both the $O(3)$ sigma model and Fendley's integrable sigma models, and we consider in
detail the effect of the theta angle on their non-perturbative corrections, by using the methods and ideas of \cite{mmr-an}.
We find that these corrections might change their sign as we move
from the theory with $\vartheta=0$ to the theory $\vartheta=\pi$. Non-perturbative corrections due to instantons
change sign according to the value of their topological charge, as expected. In particular, the leading, real exponentially small correction of the $O(3)$ sigma model changes sign as we turn on $\vartheta=\pi$, in agreement with the instanton interpretation put forward in \cite{mmr-an,bbh}. There are corrections
which are due to renormalon effects, and their imaginary part can be detected by a resurgent analysis of the perturbative series. As expected, and required by consistency, this imaginary part remains unchanged as we turn on the theta angle.
However, we find that the real part of renormalon corrections {\it does} change sign and is therefore
sensitive to the theta angle.

Our results have also some bearing on long-standing issues concerning the $O(3)$ sigma model. It has been known
for some time that the $\vartheta$ dependence of the vacuum energy and the conventional
topological susceptibility are not good observables in this model. This is best established by lattice
calculations that show that these quantities diverge in the continuum limit (see e.g.
\cite{bnpw, bietenholz} for a summary and references). An evaluation directly in the continuum with instanton
methods runs into similar problems: the integration over the instanton size is
afflicted with the usual IR divergence and an additional logarithmic UV
divergence \cite{jevicki, forster,bl,ffs} (the latter is sometimes regarded as a counterpart of the
divergence found in the lattice). Our results indicate that, in the presence of an
$h$ field, the vacuum energy as a function of $\vartheta$ is well defined (after subtracting the result for $h=0$).
The results of this paper and \cite{mmr-an} lead to a very precise prediction for the weak-coupling limit of the
resulting $h$-dependent topological susceptibility: it is given by the one-instanton effect
calculated in \cite{mmr-an}, see (\ref{chith}) for
the explicit result. It might be possible to test this prediction with a lattice calculation.

This paper is organized as follows. In section \ref{rba-sec} we consider the $O(3)$ sigma model at $\vartheta=0$,
and we review and extend results of \cite{mmr-an} on its non-perturbative corrections. In section \ref{O3-pi} we consider the
$O(3)$ sigma model at $\vartheta=\pi$, and we analyze the effect of the theta angle on non-perturbative effects. We check our analytic results with
a detailed numerical study of the Bethe ansatz equations for the massless theory at $\vartheta=\pi$. In section \ref{fendley-sec} we extend these
results to Fendley's integrable sigma models. In \ref{conclu} we list some conclusions and prospects for future research. The first appendix
gives some of the details of the computations in section \ref{rba-sec}, while the second one explains our
numerical study of the integral equation when $\vartheta=\pi$.

\sectiono{The $O(3)$ sigma model at \texorpdfstring{$\vartheta=0$}{theta = 0}}
\label{rba-sec}
\subsection{The free energy}

The observable that we will focus on in this paper is the free energy of integrable sigma models in the
presence of a magnetic field or chemical potential $h$
. The main advantage of this free energy is that,
on the one hand, it can be computed in a simple way by using the Bethe ansatz \cite{pw, wiegmann2},
and on the other hand it can be studied in perturbation theory in the regime $h \gg \Lambda$,
where $\Lambda$ is the dynamically generated scale (in this paper we will always use the
$\overline{\text{MS}}$ scheme). For this reason,
this quantity has been studied extensively in the past, and famously led to a determination of the ratio
$m/\Lambda$, first in the $O(n)$ sigma model \cite{hmn,hn}, and subsequently in
many other models \cite{pcf,fnw1,fnw2,eh-ssm,eh-scpn, hollo-lie}
(see \cite{eh-review} for a review of these results).

The free energy is defined and calculated as follows. Let $\mathsf{H}$ be the Hamiltonian of the
model, and let $\mathsf{Q}$ be the charge associated to a global conserved
current. Let $h$ be an external field coupled to $\mathsf{Q}$, which can be regarded as a chemical potential.
Let us consider the ensemble defined by the operator
\begin{equation}
\label{HQ}
\mathsf{H}- h \mathsf{Q},
\end{equation}
as well as the corresponding free energy per unit volume
\begin{equation}
\label{free-en}
F(h) =-\lim_{V, \beta \rightarrow \infty} {1\over V\beta } \log {\rm Tr}\,  {\rm e}^{-\beta (\mathsf{H}-h \mathsf{Q})},
\end{equation}
where $V$ is the volume of space and $\beta$ is the total length of Euclidean time. The quantity
\begin{equation}
\label{cfh}
 {\cal F} (h)= F(h)- F(0)
\end{equation}
can be computed by using the exact $S$ matrix and the Bethe ansatz \cite{pw,wiegmann2}. The solution is encoded in the following integral equation for a Fermi density $\epsilon(\theta)$:
\begin{equation}
\label{eps-ie}
\epsilon(\theta)-\int_{-B}^B  K(\theta-\theta') \epsilon (\theta') {\rm d} \theta' = h-m \cosh(\theta),\quad \theta \in[-B,B].
\end{equation}
Here, $m$ is the mass of the charged particles, and with a clever choice of $\mathsf{Q}$, it is directly related to the mass gap of the theory.
The kernel of the integral equation is given by
\begin{equation}
K(\theta)={1\over 2 \pi {\rm i}} {{\rm d} \over {\rm d}\theta} \log S(\theta),
\end{equation}
where $S(\theta)$ is the appropriate $S$-matrix for the scattering of the  particles charged under $\mathsf{Q}$.
The endpoints $\pm B$ are fixed by the boundary condition
\begin{equation}
\label{boundary}
\epsilon(\pm B)=0.
\end{equation}
The free energy is then given by
\begin{equation}
\label{fh-eps}
{\cal F}(h)=-{m \over 2 \pi} \int_{-B}^B \epsilon(\theta) \cosh(\theta) {\rm d} \theta.
\end{equation}
It will also be convenient to use a ``canonical" formalism and introduce the density of particles $\rho$ and energy density $e$ through
a Legendre transform of ${\cal F}(h)$,
\begin{equation}
\label{legendre}
\begin{aligned}
\rho&=-{\cal F}'(h), \\
e(\rho)-\rho h&=\mathcal{F}(h).
\end{aligned}
\end{equation}
Alternatively,
these two quantities can be obtained from the density of Bethe roots $\chi(\theta)$. This
density is supported on an interval $[-B,B]$ and satisfies the integral equation
\begin{equation}
\label{chi-ie}
\chi(\theta)-\int_{-B}^B K(\theta-\theta') \chi (\theta') {\rm d} \theta' = m \cosh \theta,
\end{equation}
where the kernel is the same one appearing in (\ref{eps-ie}). Then, we have
\begin{equation}
\label{erho}
e={m \over 2 \pi} \int_{-B}^B  \chi(\theta)  \cosh (\theta) {\rm d} \theta, \qquad \rho={1\over 2 \pi} \int_{-B}^B \chi(\theta) {\rm d} \theta.
\end{equation}
In the canonical formulation, $B$ is fixed by the value of $\rho$, and this leads implicitly to a function $e(\rho)$.

A powerful approach to solve the integral equation (\ref{eps-ie}) is the Wiener--Hopf method, as explained in \cite{jap-w,hmn,fnw1}.
We consider the Fourier transform of the kernel,
\begin{equation}
\tilde K(\omega)= \int_{\mathbb R} {\rm e}^{{\rm i} \omega \theta} K(\theta) {\rm d} \theta,
\end{equation}
and its Wiener--Hopf factorization
\begin{equation}
\label{wh-fact}
1- \tilde K (\omega)= {1\over G_+(\omega) G_-(\omega)},
\end{equation}
where $G_\pm (\omega)$ is analytic in the upper (respectively, lower) complex half plane.
Since $K(\theta)$ is an even function, we have $G_-(\omega)= G_+(-\omega)$. The result of the Wiener--Hopf analysis is the following.
We define
 %
 \begin{equation}
 \epsilon_\pm(\omega)= {\rm e}^{\pm{\rm i} B \omega} \tilde \epsilon (\omega),
 \end{equation}
 %
where $\tilde \epsilon (\omega)$ is the Fourier transform of
$\epsilon(\theta)$, after extending the latter by zero outside the
interval $[-B, B]$. Consider also the function
\begin{equation}
g_+(\omega)=  {\rm i} h \frac{ 1 - {\rm e}^{2{\rm i} B \omega} }{\omega}+
{{\rm i} m {\rm e}^B \over 2} \left( {{\rm e}^{2 {\rm i} B \omega} \over \omega-{\rm i}} -  {1\over \omega + {\rm i}}\right),
\end{equation}
which is obtained from the Fourier transform of $h - m \cosh(\theta)$, after extending it conveniently outside of $[-B,B]$, as explained in \cite{mmr-an}. Then, one finds that (\ref{eps-ie}) is equivalent to the following equation for
a function $Q(\omega)$,
\begin{equation}
\label{Q-eq}
Q(\omega)-{1\over 2 \pi {\rm i}} \int_{\mathbb R} {{\rm e}^{2 {\rm i} B\omega'} \sigma (\omega') Q(\omega') \over \omega+ \omega'+ {\rm i} 0} {\rm d} \omega'= {1\over 2 \pi {\rm i}} \int_{\mathbb R} {G_-(\omega') g_+(\omega') \over \omega+ \omega'+ {\rm i} 0} {\rm d} \omega',
\end{equation}
where
\begin{equation}
\label{sigmaeq}
\sigma(\omega) = {G_-(\omega) \over G_+(\omega)}.
\end{equation}
The solution $Q(\omega)$ to (\ref{Q-eq}) determines $\epsilon_+(\omega)$ through
\begin{equation}
\label{epsom}
{\epsilon_+(\omega) \over G_+(\omega)} = {1\over 2 \pi {\rm i}} \int_{\mathbb R} {G_-(\omega') g_+(\omega') \over  \omega'-\omega- {\rm i} 0}{\rm d} \omega'+
{1\over 2 \pi {\rm i}} \int_{\mathbb R} {{\rm e}^{2 {\rm i} B\omega'} \sigma (\omega')Q(\omega') \over\omega'- \omega- {\rm i} 0} {\rm d} \omega'.
\end{equation}
The relationship between $h$, $m$ and $B$ is determined by the boundary condition \eqref{boundary},
which in Fourier space takes the form
\begin{equation}
\label{bc-eps}
\lim_{\kappa \to +\infty} \kappa \epsilon_+({\rm i} \kappa)=0.
\end{equation}
We then find from (\ref{fh-eps}) that
\begin{equation}
\label{fh-eps-2}
{\cal F}(h)= -{1\over 2 \pi} m {\rm e}^B \epsilon_+({\rm i}).
\end{equation}
These equations were analyzed in \cite{hmn,hn} and subsequent works to obtain the perturbative expansion of ${\cal F}(h)$ in the regime $h \gg m$.
In the case of ``bosonic models'', which includes coset sigma models and their supersymmetric extensions, one can obtain a universal formula for the
one-loop free energy. In these bosonic models, the Wiener--Hopf factor $G_+({\rm i} \xi)$ has the following expansion around $\xi=0$,
\begin{equation}
\label{Gxi-exp}
G_+({\rm i} \xi)=  {k \over {\sqrt{\xi}}} \left( 1- a \xi \log \xi-b \xi+ {\cal O}(\xi^2) \right),
\end{equation}
and one finds \cite{pcf},
\begin{multline}
\label{ftba}
\mathcal{F}(h) = -{k^2 h^2 \over 4} \Biggl\{ \log\left( {h \over m} \right)+\left( a+{1\over 2} \right) \log \log\left( {h \over m} \right)+ \log \left( {{k\sqrt{2 \pi}}  \over G_+({\rm i})} \right)\\
+ a \big( \gamma_E -1+ \log(8)\big) - b -1 + \mathcal{O}\big(\log^{-1/2}(h/m)\big) \Biggr\}.
\end{multline}
A fully analytic derivation of this formula was presented in \cite{mmr-an}.

In the case of the $O(3)$ sigma model, we make the choice of charges in
\cite{wiegmann2, hmn,hn}. The resulting Euclidean Lagrangian, including the $h$ term, is
given by
\begin{equation}
\label{hLag}
{\cal L}_h = \frac{1}{2 g_0^2}\bigg\{ \partial_\mu \boldsymbol{S}\cdot \partial^\mu\boldsymbol{S}+2{\rm i} h
(S_1\partial_0 S_2-S_2 \partial_0 S_1)+h^2\left(S_3^2-1\right)\bigg\},
\end{equation}
where $\boldsymbol{S}=(S_1, S_2, S_3)$ is the quantum field satisfying the constraint $\boldsymbol{S}^2=1$,
and $g_0^2$ is the bare coupling constant. This model is
asymptotically free \cite{polyakov}. In the convention in which the beta function is given by
\begin{equation}
\label{betaf}
\beta (g)= \mu {{\rm d} g^2 \over {\rm d} \mu} = -\beta_0 g^4 - \beta_1 g^6- \cdots,
\end{equation}
one has (see e.g. \cite{bzj})
\begin{equation}
\beta_0= {1 \over 2 \pi}, \qquad \beta_1= {1\over 4 \pi^2}.
\end{equation}
In the $O(3)$ sigma model, the integral equation (\ref{eps-ie}) has the kernel
\begin{equation}
\label{kktilde}
K(\theta)={1\over \pi^2+ \theta^2}, \qquad \tilde K (\omega)= {\rm e}^{-\pi |\omega|},
\end{equation}
and its Wiener--Hopf decomposition is given by
\begin{equation}
G_+(\omega) = \frac{{\rm e}^{{\rm i}\omega\bigl[-\frac{1}{2}+\frac{\log(2)}{2}\bigr] + \frac{1}{2}{\rm i} \omega \log(-{\rm i}\omega)}}{\sqrt{-{\rm i} \omega}} \frac{\Gamma (1-{\rm i}\omega) }{\Gamma \left(\tfrac{1}{2}-\tfrac{1}{2}{\rm i}\omega\right)}.
\label{eq_G+_O3}
\end{equation}
The mass $m$ appearing in (\ref{eps-ie}) is related to the dynamically generated scale by \cite{hmn}
\begin{equation}
\label{mass-gap}
{m \over \Lambda}= {8 \over {\rm e}},
\end{equation}
and we note that
\begin{equation}
\Lambda \approx \mu \, {\rm e}^{-2 \pi/g^2(\mu)},
\end{equation}
where $g^2(\mu)$ is the renormalized coupling at the scale $\mu$.

\subsection{Non-perturbative corrections}
\label{sec_nonpert_theta_0}
Since ${\cal F}(h)$ is an observable determined by a simple integral equation, it provides an ideal testing
ground to determine both the perturbative expansion and its non-perturbative corrections. This
was first explored in the case of the principal chiral field, and in the large $N$ limit, in \cite{fkw1,fkw2}. One
possible strategy to calculate exponentially small corrections is to use the theory of resurgence combined
with perturbation theory at large orders. Volin \cite{volin, volin-thesis} developed a systematic procedure to extract long perturbative
series for ${\cal F}(h)$ directly from the Bethe ansatz equations, and this was used in
\cite{volin, volin-thesis, mr-ren,abbh1,abbh2, mmr, dpmss, bbv} to study numerically and analytically the presence of renormalon effects.

In this paper we will study non-perturbative effects in various integrable sigma models by using the method
developed in \cite{mmr-an}, which provides exact, analytic results for
the exponentially small corrections to the free energy, at finite $N$. The first non-perturbative correction in the
$O(3)$ sigma model was obtained in \cite{mmr-an} at the very first orders in the coupling constant. The calculation in \cite{mmr-an}
was verified and extended to NNLO in \cite{bbh}.

\begin{figure}[!ht]
\leavevmode
\begin{center}
\includegraphics[height=5.8cm]{figures/contouram-def.pdf}
\end{center}
\caption{The Hankel contour ${\cal C}$ can be deformed into an integral along the discontinuity of $\sigma(\omega)$,
denoted by the dashed line, plus a sum over residues. However, due to
the branch cut along the imaginary axis, this can be done in two different ways, which leads to two different integrations along the discontinuity, corresponding to the contours ${\cal C}_\pm$. The residues of the poles $\sigma_n^\pm$ will also depend on this choice.}
\label{contourdef-fig}
\end{figure}

In order to incorporate non-perturbative corrections in the Wiener--Hopf method, we deform the integration contour along the real axis in (\ref{Q-eq}) into a Hankel contour ${\cal C}$ around the positive imaginary axis. This contour is made of two rays,
one of them to the left of the imaginary axis,
and the other one to the right. If $\sigma(\omega)$ had only poles, the contour integral
could simply be evaluated by residues, as in \cite{zamo-mass} and the deformed sigma models of \cite{foz,fendley}.
However, we have to take into account the branch cut along the imaginary axis. To do this, we move the branch cut away from the imaginary axis by a small angle $\delta$. Then, as seen in \figref{contourdef-fig}, the discontinuity and the poles
become disentangled, and the integral along the path ${\cal C}$ can be separated into an integral along the discontinuity with angle $\delta$,  and a sum over the residues. The resulting tilted paths, corresponding to $\delta> 0$ (respectively, $\delta<0$) will be denoted by
${\cal C}_\pm$. As emphasized in \cite{mmr-an}, the value of the residues is sensitive to the sign of $\delta$, that is, to the branch choice of $\sigma(\omega)$. This leads to the renormalon ambiguity first discussed by David \cite{david2,david3}.
We will denote the poles of $\sigma({\rm i} \xi)$ along the positive imaginary axis by $\xi_n$. In the $O(3)$ model they can be
read from (\ref{eq_G+_O3}) and they are given by
\begin{equation}
\xi_n = n, \qquad n \in {\mathbb N}.
\end{equation}
We will also denote by $\sigma^\pm_n$ the residue of $\sigma(\omega)$ at ${\rm i}\xi_n\pm 0$, which reads
\begin{equation}
{\rm i} \sigma^\pm_{2n} = \pm  {\rm i}  \left(\frac{n}{{\rm e}}\right)^{2 n}\frac{ 2 n }{(n!)^2},\qquad {\rm i} \sigma^\pm_{2n+1} = 0.
\end{equation}
We note that the residues vanish for odd $n$, so we have removable singularities at these points\footnote{One could label only the even residues and neglect the odd ones (${\rm i} \sigma^\pm_{2n}={\rm i} \sigma'^\pm_{n}$, with $\xi'_n=2n$). However, as we will see later, the final trans-series is organized in powers of ${\rm e}^{-2nB}$ due to instanton effects. Thus, in the current numbering, $n$ tracks the exponential orders of the final trans-series.}. One finds,
\begin{equation}
\label{Q-eq-np}
Q({\rm i} \xi)-{1\over2 \pi {\rm i} } \int_{{\cal C}_\pm} { {\rm e}^{-2 B \xi'} \delta \sigma ({\rm i} \xi') Q({\rm i} \xi') \over \xi + \xi'} {\rm d} \xi'+ \sum_{n \ge 1} {{\rm e}^{-2B \xi_n} {\rm i}\sigma^\pm_n  Q_n  \over \xi + \xi_n} ={1\over 2 \pi {\rm i}} \int_{\mathbb R} {G_-(\omega') g_+(\omega' ) \over {\rm i}\xi+ \omega'+{\rm i} 0}{\rm d}\omega'.
\end{equation}
In this equation, $Q_n \equiv Q({\rm i} \xi_n)$, and
\begin{equation}
\delta \sigma ({\rm i} \xi) = \sigma(\xi({\rm i}+0)) - \sigma(\xi({\rm i}-0))
\end{equation}
is the discontinuity of $\sigma(\omega)$ across the positive
imaginary axis.
An expression similar to (\ref{Q-eq-np}), including exponentially small corrections, can be obtained from (\ref{epsom}):
\begin{equation}
\label{eps-disc}
\frac{\epsilon_+({\rm i}\kappa)}{G_+({\rm i}\kappa)} = {1\over 2 \pi {\rm i}} \int_{\mathbb R} {G_- (\omega') g_+ (\omega')
 \over \omega'-{\rm i}\kappa-{\rm i} 0}  {\rm d} \omega'+
{1\over 2 \pi {\rm i} } \int_{{\cal C}_\pm}  {{\rm e}^{-2 B \kappa'} \delta \sigma ({\rm i} \kappa') Q ({\rm i} \kappa') \over  \kappa'-\kappa} {\rm d} \kappa'
- \sum_{n \ge 1}  {{\rm e}^{-2B \xi_n} {\rm i} \sigma^\pm_n Q_n \over \xi_n-\kappa}.
\end{equation}
Note that the driving term also contains exponentially small corrections.
The free energy follows from (\ref{fh-eps-2}).

In order to study the effects of the theta angle in the $O(3)$ sigma model, it is
interesting to calculate explicitly the non-perturbative corrections up to second order in the exponentially small scale
\begin{equation}
{\rm e}^{-2 B} \approx {\rm e}^{-4 \pi/g^2(\mu)}.
\end{equation}
This calculation is very similar to what we did in \cite{mmr-an}, and we provide some of the details in Appendix \ref{app-npO3}.
Here we present the results for the free energy. One finds,
\begin{equation}
\mathcal{F}(h) = -\frac{k^2 h^2}{4} \biggl[ \mathcal{\tilde F}_{(0)}(h) +  \mathcal{\tilde F}_{(1)}(h){\rm e}^{-2B}
+ \mathcal{\tilde F}_{(2)}(h){\rm e}^{-4B} + \mathcal{O}\big({\rm e}^{-6B}\big)\biggr] \mp {\rm i} \frac{m^2}{16},
\label{eq_free_energy_exp4B}
\end{equation}
where
\begin{align}
\mathcal{\tilde F}_{(0)}(h) &= B + \frac{\log(B)- 2 + 4\log(2)}{2} + \mathcal{O}\big(B^{-1}\big),\\
\mathcal{\tilde F}_{(1)}(h) &= -\frac{4 B^2}{{\rm e}}+\frac{-2 \log(B) - 3 + 2 \gamma_E -6 \log (2)}{{\rm e}}B + \mathcal{O}\big(B^0\big),\\
\mathcal{\tilde F}_{(2)}(h) &= \frac{8 B^2}{{\rm e}^2} + \frac{4 \log(B) + 2 - 4 \gamma_E + 12\log(2) \mp 2 {\rm i} }{{\rm e}^2}B + \mathcal{O}\big(B^0\big).
\end{align}
The imaginary ambiguity in the ${\rm e}^{-4B}$ correction originates from the residue $\sigma_2^\pm$. The
expression (\ref{eq_free_energy_exp4B}) is an example of a trans-series (see \cite{mmlargen,ss,abs} for an overview),
and it involves the small parameters $1/B$, $\log(B)/B$, and ${\rm e}^{-2B}$.

As pointed out in \cite{bbbkp}, it is useful to express the free energy in terms of the coupling $\tilde\alpha$, defined by
\begin{equation}
\frac{1}{\tilde\alpha} + \log(\tilde\alpha) = \log\left( \frac{h}{\Lambda} \right).
\label{eq_alphatilde_def}
\end{equation}
It is easy to see that, at leading order, one has
\begin{equation}
\tilde\alpha \sim  \beta_0 g^2(h),
\end{equation}
where $g^2(h)$ is the running coupling constant
at the scale $h$. One uses the boundary condition (\ref{boco}) and the mass gap (\ref{mass-gap})
to obtain $B$ as a trans-series in $\tilde \alpha$.
Then one can write \eqref{eq_free_energy_exp4B} as
\begin{equation}
\label{fh-ts}
\mathcal{F}(h) = -\frac{ h^2}{4 \pi} \biggl[ \mathcal{F}_{(0)}(h) +  \mathcal{F}_{(1)}(h){\rm e}^{-2/\tilde\alpha} + \mathcal{F}_{(2)}(h){\rm e}^{-4/\tilde\alpha} + \mathcal{O}\big({\rm e}^{-6/\tilde\alpha}\big)\biggr] \mp {\rm i} \frac{m^2}{16},
\end{equation}
where
\begin{equation}
\label{f1h}
\begin{aligned}
\mathcal{F}_{(0)}(h) &= \frac{1}{\tilde\alpha} - \frac{1}{2} + \mathcal{O}\big(\tilde\alpha\big),\\
\mathcal{F}_{(1)}(h) &= -\frac{64}{{\rm e}^2 \tilde\alpha^3} + \frac{32 (-\log(\tilde\alpha) -3 + \gamma_E +5\log (2))}{{\rm e}^2 \tilde\alpha^2}+O\big(\tilde\alpha^{-1}\big),\\
\mathcal{F}_{(2)}(h) &= \frac{512(1 \mp {\rm i})}{{\rm e}^4 \tilde\alpha^3} + \mathcal{O}\big( \tilde\alpha^{-2} \big).
\end{aligned}
\end{equation}
Once again, the imaginary ambiguity in the ${\rm e}^{-4/\tilde\alpha}$ correction originates from $\sigma_2^\pm$.

Finally, we write the result in the canonical formulation, i.e. in terms of the density of particles $\rho$
and energy density $e$ introduced in (\ref{legendre}). As noted in \cite{bbbkp},
the appropriate quantity to study is the quotient $e/(\pi\rho^2)$.
In addition, we will express the result in terms of the coupling $\alpha$, defined by
\begin{equation}
\frac{1}{\alpha} = \log\left( \frac{\rho}{2\beta_0 \Lambda} \right).
\label{eq_alpha_def}
\end{equation}
One finds that $\tilde \alpha$ is given by a trans-series in $\alpha$, with $\tilde \alpha \sim \alpha$ at the very leading order.
Our final result for the normalized energy density, up to second order in the exponential corrections, is
\begin{multline}
\frac{e}{\pi\rho^2} = \alpha+\frac{\alpha^2}{2} + \mathcal{O}\big(\alpha^3\big)
+ \frac{32}{{\rm e}^2}\left( \frac{2}{\alpha } + \log(\alpha ) +3  - \gamma_E - 5\log(2) \mp {\rm i} \frac{\pi}{2}
+\frac{\alpha}{2}
+\mathcal{O}\big(\alpha^2\big) \right)  {\rm e}^{-2/\alpha}\\
+ \frac{512}{{\rm e}^4}\left( \frac{1 \pm {\rm i}}{\alpha} + \mathcal{O}\big(\alpha^0\big) \right) {\rm e}^{-4/\alpha} + \mathcal{O}\big( {\rm e}^{-6/\alpha} \big).
\label{eq_normalized_energy}
\end{multline}
The first ambiguous imaginary term
arises from the isolated $m^2$ term in \eqref{eq_free_energy_exp4B}, while  the term of order $\alpha$ in the correction proportional to ${\rm e}^{-2/\alpha}$ was calculated in \cite{bbh}\footnote{As this paper was being typed, \cite{bbhv22} appeared, which calculates this series up to order $\alpha^3$, and estimates numerically the coefficients up to order $\alpha^5$.}. The expression (\ref{eq_normalized_energy}), up to and including the first exponential correction, was
tested in detail in \cite{mmr-an,bbh} against a numerical calculation based on the canonical formalism of (\ref{chi-ie}), (\ref{erho}).

\section{The $O(3)$ sigma model at \texorpdfstring{$\vartheta =\pi$}{theta = pi}}
\label{O3-pi}
As we mentioned in the introduction, the $O(3)$ sigma model admits a theta angle. The topological density is of the form
\begin{equation}
q(x)= {1\over 8 \pi} \epsilon^{abc} \epsilon_{\mu \nu} S_a \partial_\mu S_b \partial_\nu S_c,
\end{equation}
where $\epsilon^{abc}$, $\epsilon_{\mu \nu}$ are totally antisymmetric symbols for the internal and spacetime indices,
respectively. The topological charge is then given by
\begin{equation}
\label{top-charge}
{\cal Q}=\int q(x) {\rm d}^2 x,
\end{equation}
 and one can then add to the Euclidean action the term
\begin{equation}
{\rm i} \vartheta  {\cal Q}
\end{equation}
The $\vartheta$ angle has important effects on the theory. Haldane famously conjectured \cite{haldane} that,
when $\vartheta= \pi$, the mass gap vanishes and the sigma model flows to a non-trivial
conformal field theory, the $SU(2)_1$ Wess--Zumino--Witten model
(see e.g. \cite{affleck-review} for a review of these results). In particular, we can consider the free energy (\ref{cfh}) in the presence of the
theta angle,
\begin{equation}
\label{cfh-theta}
{\cal F}(h, \vartheta) = F(h, \vartheta)- F(0, \vartheta),
\end{equation}
and we can ask how it changes as we vary $\vartheta$. Unfortunately, it is very difficult to answer this question
analytically for arbitrary values of $\vartheta$. However, when $\vartheta=\pi$, the $O(3)$ non-linear sigma model is expected to be integrable,
and its massless scattering theory
has been proposed in \cite{zz-theta}, building on \cite{fz}. This
makes it possible to calculate ${\cal F}(h, \pi)$. In this section we will present some detailed
results for this quantity, with emphasis on its non-perturbative corrections.

\subsection{The free energy}

The general setting to calculate the free energy in a massless scattering theory has been developed in
\cite{foz, fsz1, fsz2} (see \cite{fs-rev} for a review). We consider again a conserved charge, coupled to a magnetic field or
chemical potential $H$. The ground state of the massless theory
is described by two rapidity distributions $\epsilon_1$ and $\epsilon_2$, which
characterize right and left moving particles, respectively. The domain of these
distributions is only bounded on one-side, since massless particle can
have arbitrarily low momentum. The boundary condition is then
\begin{equation}
\epsilon_1(B) = \epsilon_2(-B) = 0.
\label{bceq}
\end{equation}
The Bethe ansatz leads to a pair of integral equations for the rapidity distributions $\epsilon_{1,2}$, which have the form
\begin{equation}
\begin{aligned}
\epsilon_1(\theta) - \int_{-\infty}^B \varphi_1(\theta-\theta')\epsilon_1(\theta'){\rm d} \theta' - \int^{\infty}_{-B} \varphi_2(\theta-\theta')\epsilon_2(\theta'){\rm d} \theta' &= H - \frac{M {\rm e}^\theta}{2}, \quad &\theta&<B,\\
\epsilon_2(\theta)  - \int_{-\infty}^B \varphi_2(\theta-\theta')\epsilon_1(\theta'){\rm d} \theta' - \int^{\infty}_{-B} \varphi_1(\theta-\theta')\epsilon_2(\theta'){\rm d} \theta' &= H - \frac{M {\rm e}^{-\theta}}{2}, \quad &\theta&> -B .
\end{aligned}
\label{orig_elr_inteq}
\end{equation}
In the cases we will consider, the two functions are related by parity
$\epsilon_1(\theta) = \epsilon_2 (-\theta)$. The kernels
$\varphi_1$ and $\varphi_2$ are extracted respectively from the left-left (or right-right)
massless $S$-matrix and the left-right (or right-left) massless $S$-matrix
corresponding to the particles charged under $\mathsf{Q}$. In these equations, $M$ is
a mass scale related to the mass gap $m$ of the theory at
$\vartheta=0$. With the above normalization,
we have $M=m$. In the case of the massless theory for the $O(3)$ non-linear sigma model at $\vartheta=\pi$,
the chemical potential $H$ is related to the chemical potential $h$ at $\vartheta =0$ simply by
\begin{equation}
\label{Hhr}
H={h \over r},
\end{equation}
where as shown e.g. in \cite{fs-rev} the constant $r$ can be obtained from the kernels in (\ref{orig_elr_inteq}),
see (\ref{hH}) below for an explicit formula. The free energy is then given by
\begin{equation}
{\cal F}(h, \pi)= -{M \over 2 \pi} \int_{-\infty}^B {\rm e}^\theta \epsilon_1(\theta) {\rm d} \theta.
\label{free-en-massless}
\end{equation}

As in the massive phase at $\vartheta=0$, it is useful to consider the density of particles and energy introduced in (\ref{legendre}). These can be obtained,
as in (\ref{erho}), (\ref{chi-ie}), from ``dual"
integral equations for Bethe root densities $\chi_1(\theta)$, $\chi_2(\theta)$,
\begin{equation}
\label{merTBA}
\begin{aligned}
\chi_1 (\theta)- \int_{-\infty}^B \varphi_1 (\theta-\theta') \chi_1(\theta'){\rm d} \theta'-  \int_{-B}^\infty \varphi_2 (\theta-\theta') \chi_2(\theta'){\rm d} \theta' &= M {\rm e}^\theta, \\
\chi_2 (\theta)- \int_{-B}^\infty \varphi_1 (\theta-\theta') \chi_2(\theta'){\rm d} \theta'-  \int_{-\infty}^B \varphi_2 (\theta-\theta') \chi_1(\theta'){\rm d} \theta' &= M  {\rm e}^{-\theta}.
\end{aligned}
\end{equation}
We can bring these two integral equations into a single equation for $\chi_1(\theta)$, using the parity relation $\chi_1(-\theta) = \chi_2(\theta)$ and the symmetry of the kernels  $\varphi_i(\theta) = \varphi_i(-\theta)$:
\begin{equation}
\chi_1(\theta) - \int_{-\infty}^B \big( \varphi_1(\theta-\theta') + \varphi_2(\theta + \theta') \big) \chi_1(\theta') {\rm d} \theta' =  M{\rm e}^\theta.
\label{merTBA2}
\end{equation}
Then the density of particles and the energy density are given by
\begin{equation}
e={ M \over 4 \pi} \int_{-\infty}^B {\rm e}^\theta \chi_1 (\theta) {\rm d} \theta, \qquad \rho={ 1 \over 2 \pi r } \int_{-\infty}^B \chi_1 (\theta) {\rm d} \theta.
\label{erho_pi}
\end{equation}

\subsection{Non-perturbative corrections}
\label{sec_O3pi}

We can now try to extract both the perturbative part and the non-perturbative corrections from ${\cal F}(h, \pi)$ in the $O(3)$ sigma model
with $\vartheta = \pi$ and compare them to what we found for $\vartheta=0$. We expect that the perturbative part remains unchanged, and
we would like to understand what happens to the non-perturbative effects. These issues were already addressed
in e.g. \cite{foz,fendley}. However, in their calculation of the free energy, \cite{foz,fendley} considered {\it deformations} of the conventional sigma
models (including the so-called ``sausage" model deformation of the $O(3)$ sigma model). These deformations
have an interest of their own, and they are relatively easier to study: after the deformation, the calculation of perturbative and
non-perturbative corrections becomes straightforward, as in the Sine-Gordon model studied in \cite{zamo-mass}.
This is due to the fact that the analogue of the function (\ref{sigmaeq})
has only pole singularities, and no branch cuts. However, the deformed
theories are not necessarily good guides for the undeformed theories. For example, in order to obtain the
corresponding results for ${\cal F}(h, \pi)$ in the undeformed theory, one should perform a highly non-trivial resummation.
In addition, due to the absence of branch cuts, the deformed theories do not seem to have renormalons, hence an important aspect of the physics
of the original models is lost.

In spite of these differences, some aspects of the analysis of the deformed theories also
apply to the undeformed theories considered in this paper.
As in \cite{foz,fendley,fsz2,fs-rev}, the first step is a Wiener-Hopf analysis
of the integral equations (\ref{orig_elr_inteq}), similar to what is done when $\vartheta=0$. Let us recall that one
of the basic ingredients in this analysis is the
additive decomposition into a function analytic in the complex upper half plane $\mathbb{H}_+$
(which carries a $+$ subscript) and a function analytic in the complex lower half plane $\mathbb{H}_-$
(which carries a $-$ subscript). For a generic function $\psi(\omega)$ defined in the real line, this can be obtained from
\begin{equation}
[\psi(\omega)]_\pm =
 \pm \frac{1}{2\pi{\rm i}}\int_{\mathbb R} \frac{\psi(\omega')}{\omega'-(\omega\pm{\rm i} 0)}{\rm d}\omega',\quad \omega\in \mathbb{H}_\pm.
\label{WH_decompostion}
\end{equation}
After analytically extending each factor to the entire plane (up to singularities), we have
\begin{equation}
\psi(\omega) = [\psi(\omega)]_+ + [\psi(\omega)]_-.
\end{equation}
A useful consequence of \eqref{WH_decompostion} is that if a function $\Psi_\pm(\omega)$ is analytic and bounded in the upper (lower) half plane, then
\begin{equation}
[\Psi_\pm]_\mp = 0, \qquad [{\rm e}^{\pm  {\rm i} a \omega}\Psi_\pm]_\mp = 0,\quad a>0.
\label{WH_cancel}
\end{equation}

The goal of this section is to rewrite the integral equations \eqref{orig_elr_inteq} to make the comparison with the $\vartheta=0$ case manifest. Thus we will first reduce them to individual integral equations of the form of \eqref{Q-eq-np} and \eqref{eps-disc} and then single out the perturbative contributions, the non-perturbative ambiguous formally imaginary contributions and the non-perturbative formally real unambiguous contributions. As we previously discussed, we expect the first to be identical in the $\vartheta=0$ and $\vartheta=\pi$ cases. The picture is slighly more complicate for non-perturbative corrections. Roughly speaking, imaginary ambiguous terms will remain identical, which is consistent with their connection to large order behaviour of the perturbative series, while formally real non-perturbative effects will change sign according to the rule ${\rm e}^{-2nB}\rightarrow (-1)^n {\rm e}^{-2nB}$. This description is complicated by the fact that real and imaginary exponential terms will sometimes multiply each other when constructing the free energy from  $Q({\rm i}\xi)$ and $\epsilon_+({\rm i}\kappa)$. We will give a more precise statement by the end of this section.

While we also often omit the qualification, it is important to draw attention to the fact that these terms are only \textit{formally} real or imaginary, i.e.
we can say that they are real or imaginary in the sense that they are real or imaginary terms of a formal power series. However, these divergent formal
power series are only made sense of once we apply Borel summation, which combines formally real and imaginary terms into a final real, unambiguous result.

To start the Wiener--Hopf analysis of the integral equations (\ref{orig_elr_inteq}), we have to extend their
domain of validity to the full real line. We first extend the driving terms as
\begin{equation}
g(\theta) = \begin{cases}
      H- \dfrac{M{\rm e}^\theta}{2} & \theta < B, \\[2mm]
      H - h  &  \theta>B,
\end{cases}
\label{gextend}
\end{equation}
and introduce an unknown function $Y(\theta)$ which is supported in the positive half-line. The constant $h$  in \eqref{gextend} is at the moment
arbitrary and changing it amounts simply to redefining the unknown function $Y$. Then we have
\begin{equation}
\begin{aligned}
\epsilon_1(\theta) - \int_{-\infty}^B \varphi_1(\theta-\theta')\epsilon_1(\theta'){\rm d} \theta' - \int^{\infty}_{-B} \varphi_2(\theta-\theta')\epsilon_2(\theta'){\rm d} \theta' &= g(\theta) + Y(\theta-B) , \\
\epsilon_2(\theta) - \int_{-\infty}^B \varphi_2(\theta-\theta')\epsilon_1(\theta'){\rm d} \theta' - \int^{\infty}_{-B} \varphi_1(\theta-\theta')\epsilon_2(\theta'){\rm d} \theta'&= g(-\theta)+ Y(-\theta-B) .
\end{aligned}
\label{extended_elr_inteq}
\end{equation}
To find the free energy of \eqref{free-en-massless}, we will first reduce the two above equations into a single equation for $\epsilon_1$. To this end, we take the Fourier transform of the integral equations, which gives
\begin{align}
\tilde{\epsilon}_1-\phi_1 \tilde{\epsilon}_1- \phi_2 \tilde{\epsilon}_2&=
 {\rm e}^{{\rm i}  B \omega}Y_++{\rm e}^{{\rm i}  B \omega} g^m_- + H 2\pi \delta - h\frac{{\rm i}{\rm e}^{{\rm i} B \omega}}{\omega+{\rm i} 0},\label{epsilon1_fourier}\\
 \tilde{\epsilon}_2-\phi_2 \tilde{\epsilon}_1- \phi_1 \tilde{\epsilon}_2&={\rm e}^{-{\rm i}  B \omega}Y_-+{\rm e}^{-{\rm i}  B \omega} g^m_+  + H 2\pi \delta + h\frac{{\rm i}{\rm e}^{-{\rm i} B \omega}}{\omega-{\rm i} 0},
\label{epsilon2_fourier}
\end{align}
where $\tilde\epsilon_i$, $\phi_i$ and $Y_+$ are the Fourier transforms of $\epsilon_i$, $\varphi_i$ and $Y$, respectively. The last three terms in each equation correspond to the Fourier transform of $g(\theta)$ and $g(-\theta)$, respectively, where we have conveniently isolated the Fourier transform of the term proportional to $M$ in \eqref{gextend}:
\begin{equation}
g_-^{m}(\omega) = {\rm e}^{-{\rm i} B \omega}\int_{-\infty}^B {\rm e}^{{\rm i}\omega\theta} \left( -\frac{m  {\rm e}^\theta}{2} \right) {\rm d} \theta=
\frac{m {\rm e}^B}{2} \frac{{\rm i}}{\omega - {\rm i}}.
\end{equation}
From now on, we will use that $M=m$. Finally, $Y_-(\omega) = Y_+(-\omega)$, $g_+^m(\omega) = g_-^m(-\omega)$ and $\delta$ is the Dirac delta function. We can now solve for $\tilde{\epsilon}_2$ in \eqref{epsilon2_fourier}, which gives
\begin{equation}
\tilde{\epsilon}_2 = \frac{1}{1-\phi_1}\left(\phi_2 \tilde{\epsilon}_1+{\rm e}^{-{\rm i}  B \omega} Y_- + {\rm e}^{-{\rm i}  B \omega} g^m_+ + H 2\pi \delta+ h  \frac{{\rm i}{\rm e}^{-{\rm i} B \omega}}{\omega-{\rm i} 0}\right).
\end{equation}
If we plug this expression back into \eqref{epsilon1_fourier}, we obtain
\begin{multline}
\label{yint}
\left(1-\phi_1 - \frac{\phi_2^2}{1-\phi_1}\right) \epsilon_- = Y_+ + g^m_- + {\rm e}^{-2{\rm i} B \omega}\frac{\phi_2}{1-\phi_1}\left(Y_-+g^m_+\right)\\
+ H 2\pi \delta\left(1+ \frac{\phi_2(0)}{1-\phi_1(0)}\right)
- h \left(\frac{{\rm i}}{\omega+{\rm i} 0} - \frac{\phi_2}{1-\phi_1}\frac{{\rm i}{\rm e}^{-2{\rm i} B \omega}}{\omega-{\rm i} 0}\right).
\end{multline}
where
\begin{equation}
\epsilon_-(\omega) = {\rm e}^{-{\rm i}  B \omega} \tilde{\epsilon}_1(\omega).
\end{equation}
The l.h.s. of (\ref{yint}) suggests to
consider the Wiener--Hopf factorization of the kernel
\begin{equation}
\label{eff-kernel}
1-\phi_1(\omega) - \frac{\phi_2^2(\omega)}{1-\phi_1(\omega)}
= \frac{1}{K_+(\omega)K_-(\omega)},
\end{equation}
where $K_\pm(\omega)$ is analytic in the upper (respectively, lower) complex half plane. We now use the additive decomposition of the $\delta$-function
\begin{equation}
2\pi \delta(\omega) = \frac{{\rm i}}{\omega+{\rm i} 0} - \frac{{\rm i}}{\omega-{\rm i} 0}
\label{delta_def}
\end{equation}
and we see that the choice
\begin{equation}
\label{hH}
h  =  H \left(1+\frac{\phi_2(0)}{1-\phi_1(0)}\right)
\end{equation}
simplifies the terms singular at the origin. This $h$ will be eventually identified with the chemical potential of the
theory with $\vartheta=0$, therefore (\ref{hH}) provides the relation (\ref{Hhr}) stated above. Lastly, we define
\begin{equation}
v = K_+ Y_+,
\label{vdef_final}
\end{equation}
which plays the role of the function $Q$ appearing in (\ref{Q-eq}).
Algebraic manipulations then give us
\begin{equation}
\begin{aligned}
\frac{\epsilon_-}{K_-} &= v(\omega) + {\rm e}^{-2{\rm i} B \omega}\frac{\phi_2}{1-\phi_1} \sigma(-\omega) v(-\omega) + K_+\left(g^m_- +  {\rm e}^{-2{\rm i} B \omega} \frac{\phi_2}{1-\phi_1}g^m_+\right) \\
&- \frac{{\rm i} h K_+}{\omega-{\rm i} 0}  \left(1 - {\rm e}^{-2{\rm i} B \omega} \frac{\phi_2}{1-\phi_1}\right),
\label{completeeq_genform}
\end{aligned}
\end{equation}
where
\begin{equation}
\sigma(\omega) = \frac{K_-(\omega)}{K_+(\omega)}.
\end{equation}

So far our considerations have been very general. Let us now focus on the specific example
of the $O(3)$ sigma model with $\vartheta=\pi$.
After turning on the chemical potential as in (\ref{hLag}), one finds a pair of equations of the form (\ref{orig_elr_inteq}), where the Fourier
transforms of the kernels are given by \cite{foz}
\begin{equation}
\phi_1(\omega)=\phi_2(\omega) =
{1\over 1+ {\rm e}^{\pi |\omega|}}.
\end{equation}
It turns out that
\begin{equation}
1-\phi_1(\omega) - \frac{\phi_2^2(\omega)}{1-\phi_1(\omega)} = 1- \tilde K(\omega),
\label{phitoK}
\end{equation}
where $\tilde K (\omega)$ is defined in (\ref{kktilde}). Therefore, $K_\pm (\omega)= G_\pm (\omega)$,
where $G_+(\omega)$ is given in (\ref{eq_G+_O3}). Furthermore, we find the relation
\begin{equation}
\frac{\phi_2}{1-\phi_1}=1-\frac{1}{G_+G_-}.
\label{phi2phi1_simplification}
\end{equation}
We now take the $(+)$-projection of \eqref{completeeq_genform}. Due to the simplification \eqref{phi2phi1_simplification} we can cancel many terms using the property \eqref{WH_cancel},
and we find
\begin{equation}
v = -\left[ G_+ g_- \right]_+
-\left[{\rm e}^{-2{\rm i} B \omega} \left\{\sigma v\right\}(-\omega)\right]_+-\frac{{\rm i} m {\rm e}^{-B}}{2 G_+({\rm i})}  \frac{1}{\omega+{\rm i}}.
\label{v_O3}
\end{equation}
where we have introduced
\begin{equation}
g^{(p)}_+(\omega) = {\rm i} h \frac{1-{\rm e}^{2{\rm i}  B \omega}}{\omega} - \frac{{\rm i} m {\rm e}^B}{2}\frac{1}{\omega+{\rm i}}
\end{equation}
and
\begin{equation}
g_+(\omega)  =g_-(-\omega)  = g^{(p)}_+(\omega) + \frac{{\rm i} m {\rm e}^B}{2}\frac{{\rm e}^{2{\rm i} B \omega}}{\omega-{\rm i}}.
\label{gdef}
\end{equation}
Note that, although we keep a subscript $+$ for both, $g^{(p)}_+$ is analytic in the upper half plane, but $g_+$ is not.

One can also obtain an equation for $\epsilon_-$ from \eqref{completeeq_genform}.
After taking the $(-)$-projection, we find
\begin{equation}
\frac{\epsilon_-}{G_-} =
\left[
G_+ g_-
\right]_-
+ \left[
{\rm e}^{-2{\rm i} B \omega}
\left\{\tilde\sigma v \right\}(-\omega)
\right]_- +
\frac{{\rm i} m {\rm e}^B}{2(\omega+{\rm i})}
\left(\frac{{\rm e}^{-2{\rm i}  B \omega}}{G_-(\omega)}- \frac{{\rm e}^{-2 B}}{G_+({\rm i})}
\right)
- \frac{{\rm i} h{\rm e}^{-2{\rm i} B \omega}}{\omega G_-(\omega)},
\label{eminus_O3}
\end{equation}
where we now introduce
\begin{equation}
\tilde\sigma(\omega)=\frac{G_-(\omega)}{G_+(\omega)}-\frac{1}{G_+^2(\omega)}.
\label{atilde}
\end{equation}
We remark that the extra term $1/G_+^2(\omega)$ is analytic in the upper half plane, and this will have important consequences. In particular, it will imply that the perturbative part does not change when going from $\vartheta=0$ to $\vartheta=\pi$.

In order to make contact with the analysis of \cite{mmr-an}, it is useful to write \eqref{v_O3} and \eqref{eminus_O3} in
integral form. For $v(\omega)$ we have
\begin{equation}
v({\rm i} \xi) = \frac{1}{2\pi {\rm i}}\int_{{\mathbb R}} \frac{G_-(\omega')g_+(\omega')}{{\rm i}\xi+\omega'+{\rm i} 0}{\rm d} \omega' +\frac{1}{2\pi {\rm i}}\int_{{\mathbb R}} \frac{{\rm e}^{2{\rm i} B \omega'}\sigma(\omega')v(\omega')}{{\rm i}\xi+\omega'+{\rm i} 0}  {\rm d} \omega' - \frac{m {\rm e}^{-B}}{2G_+({\rm i})}\frac{1}{\xi+1},
\label{veq}
\end{equation}
and for $\epsilon$ we get
\begin{multline}
\frac{\epsilon_+({\rm i}\kappa)}{G_+({\rm i}\kappa)}=
\frac{1}{2\pi{\rm i}}\int_{{\mathbb R} } \frac{G_-(\omega')g_+(\omega')}{\omega'-{\rm i}\kappa-{\rm i} 0} {\rm d} \omega'
+\frac{1}{2\pi{\rm i}}\int_{{\mathbb R} }\frac{{\rm e}^{2{\rm i} B\omega'} \tilde\sigma(\omega')v(\omega')}{\omega'-{\rm i}\kappa-{\rm i} 0} {\rm d} \omega' \\
-\frac{m{\rm e}^{B}}{2(\kappa-1)}\left(\frac{{\rm e}^{-2 B \kappa}}{G_+({\rm i}\kappa)}-\frac{{\rm e}^{-2 B}}{G_+({\rm i})}\right)
+ \frac{h {\rm e}^{-2 B \kappa}}{\kappa G_+({\rm i} \kappa)}.
\label{eps_eq}
\end{multline}
Equations (\ref{veq}) and (\ref{eps_eq}) are respectively the counterparts in the massless
theory of equations (\ref{Q-eq}) and (\ref{epsom}) in the massive theory. We can analyze them
as in \cite{mmr-an}, and extract the non-perturbative corrections by deforming the contour as we did for the theory at $\vartheta=0$.

Let us first consider \eqref{veq}. The two integrals contain perturbative terms due to the discontinuity, non-perturbative imaginary terms due to the poles of $G_-$ (and thus $\sigma$) while only the first integral has a non-perturbative real contribution from the pole of $g_+$ at $\omega={\rm i}$. We want to single out this last contribution. Since $g^{(p)}_+$ does not have such a pole, we expand the first integral in \eqref{veq} and
 plug in
\begin{equation}
G_-({\rm i}) = \frac{1}{\sqrt{2{\rm e}}},\qquad G_+({\rm i})= \sqrt{\frac{{\rm e}}{2}}.
\label{K+K+}
\end{equation}
We then find
\begin{multline}
v({\rm i} \xi) = \frac{1}{2\pi {\rm i}}\int_{\mathbb R} \frac{ G_-(\omega')g^{(p)}_+(\omega')}{{\rm i}\xi+\omega'+{\rm i} 0}  {\rm d} \omega'
+\frac{1}{2\pi {\rm i}}\int_{\mathbb{R}}  \frac{{\rm e}^{2{\rm i} B \omega'}\sigma(\omega')v(\omega')}{{\rm i}\xi+\omega'+{\rm i} 0} {\rm d} \omega'
\\
+\frac{m{\rm e}^B}{4\pi{\rm i}}\int_{{\cal C}_\pm}  \frac{ {\rm e}^{-2B\xi'}\delta G_-({\rm i}\xi')}{(\xi+\xi')(\xi'-1)}{\rm d} \xi'
- \frac{m{\rm e}^B}{2}\sum_{n\ge 1} \frac{{\rm e}^{-2B\xi_n} G_+({\rm i}\xi_n){\rm i}\sigma_n^\pm}{(\xi_n+\xi)(\xi_n-1)}
- \frac{ m {\rm e}^{-B}}{2\sqrt{2{\rm e}}}\frac{1}{\xi+1} \,.
\label{veq_expand}
\end{multline}
In this form it is easier to break down the several contributions. Since $g^{(p)}_+$ and $v$ are analytic in the upper half plane, we can deform the first two integrals around the positive imaginary axis, which returns a perturbative part due to the discontinuity, and non-perturbative imaginary terms due to the poles. The third term in \eqref{veq_expand} is already an integral along the discontinuity which can be expanded into a perturbative series. The terms $G_+({\rm i}\xi_n)\sigma_n^\pm$ in the sum are imaginary and ambiguous, while the last term in \eqref{veq_expand} is the sole real unambiguous non-perturbative term.

For $\vartheta=0$, a similar expansion of \eqref{Q-eq-np} gives the following analogous formula:
\begin{multline}
Q({\rm i} \xi) = \frac{1}{2\pi {\rm i}}\int_{\mathbb R} \frac{ G_-(\omega')g^{(p)}_+(\omega')}{{\rm i}\xi+\omega'+{\rm i} 0} {\rm d} \omega'
+\frac{1}{2\pi {\rm i}}\int_{\mathbb{R}} \frac{{\rm e}^{2{\rm i} B \omega'}\sigma(\omega')Q(\omega')}{{\rm i}\xi+\omega'+{\rm i} 0}  {\rm d} \omega'
\\
+\frac{m{\rm e}^B}{4\pi{\rm i}}\int_{{\cal C}_\pm}  \frac{ {\rm e}^{-2B\xi'}\delta G_-({\rm i}\xi')}{(\xi+\xi')(\xi'-1)} {\rm d} \xi'
- \frac{m{\rm e}^B}{2}\sum_{n\ge 1} \frac{{\rm e}^{-2B\xi_n} G_+({\rm i}\xi_n){\rm i}\sigma_n^\pm}{(\xi_n+\xi)(\xi_n-1)}
+ \frac{ m {\rm e}^{-B}}{2\sqrt{2{\rm e}}}\frac{1}{\xi+1} \,.
\label{Q_expand}
\end{multline}
Comparing \eqref{veq_expand} and \eqref{Q_expand} with $Q\leftrightarrow v$, we see that only
the real exponentially suppressed contribution has the opposite sign. While it is not apparent from this derivation, it is useful to keep in mind that
\begin{equation}
    h\sim m {\rm e}^{B},
\end{equation}
and thus $m {\rm e}^{-B}$ is a correction proportional to ${\rm e}^{-2B}$ once it has been expressed in terms of $h$.

Let us do a similar analysis of the equation for $\epsilon_+({\rm i}\kappa)$, \eqref{eps_eq}. This equation is used in the calculation of the free energy in two distinct ways. First, to compute the free energy proper in \eqref{free-en-massless} from the Wiener-Hopf solution. In this case, we evaluate \eqref{eps_eq} at $\kappa=1$ to obtain
\begin{equation}
{\cal F}(h, \pi) = - \frac{1}{2\pi} m{\rm e}^B \epsilon_+({\rm i}).
\end{equation}
Second, to compute the relation between $B$ and $m/h$. In this case, we impose the boundary condition \eqref{bceq}, which in Fourier space takes the form:
\begin{equation}
\lim_{\kappa\rightarrow\infty}\kappa\epsilon_+({\rm i}\kappa)=0,
\end{equation}
Thus, we have to compute the limit of \eqref{eps_eq} as $\kappa \rightarrow \infty$. These two evaluations of $\epsilon_+({\rm i}\kappa)$, at $\kappa=1$ and in the limit $\kappa\rightarrow\infty$, must be analyzed separately.

Let us start with the $\kappa=1$ case by observing that the second term in the r.h.s. of \eqref{eps_eq} can be written as
\begin{equation}
\begin{aligned}
\frac{1}{2\pi{\rm i}}\int_{\mathbb R} \frac{{\rm e}^{2{\rm i} B\omega'} \tilde\sigma(\omega')v(\omega')}{\omega'-{\rm i}-{\rm i} 0} {\rm d} \omega'  &= \frac{1}{2\pi{\rm i}}\int_{{\cal C}_\pm} \frac{{\rm e}^{-2 B\xi'} \delta\sigma({\rm i}\xi')v({\rm i}\xi')}{\xi'-1} {\rm d} \xi' - \sum_{n\ge 1} \frac{{\rm e}^{-2B\xi_n}{\rm i}\sigma_n^\pm v({\rm i}\xi_n)}{\xi_n-1} \\
&  - \frac{{\rm e}^{-2B}  }{{\rm e}}v({\rm i}),
\end{aligned}
\end{equation}
where we use that discontinuities and residues of $\tilde{\sigma}$ are identical to those of $\sigma$,
since they differ only by the addition of the function $1/G_+^2(\omega)$, which is analytic in the upper half plane, as we remarked in \eqref{atilde}. This formula contrasts with the
analogous integral in the $\vartheta=0$ calculation:
\begin{equation}
\begin{aligned}
\frac{1}{2\pi{\rm i}}\int_{\mathbb R}  \frac{{\rm e}^{2{\rm i} B\omega'} \sigma(\omega')Q(\omega')}{\omega'-{\rm i}-{\rm i} 0} {\rm d} \omega' &= \frac{1}{2\pi{\rm i}}\int_{{\cal C}_\pm} \frac{{\rm e}^{-2 B\xi'} \delta\sigma({\rm i}\xi')Q({\rm i}\xi')}{\xi'-1} {\rm d} \xi' - \sum_{n\ge 1} \frac{{\rm e}^{-2B\xi_n}{\rm i}\sigma_n^\pm Q({\rm i}\xi_n)}{\xi_n-1} \\
& + \frac{{\rm e}^{-2B}  }{{\rm e}}Q({\rm i}).
\end{aligned}
\end{equation}
Thus for this term we have once again that the perturbative and imaginary non-perturbative terms are identical, but the real non-perturbative term changes sign.

Furthermore, in the first integral of \eqref{eps_eq}, we also have real exponentially suppressed terms, which come from the residue of the double pole at $\omega = {\rm i}$, and they are identical to the $\vartheta=0$ case. Explicitly, we find
\begin{multline}
\label{G-g+_computation}
\frac{1}{2\pi{\rm i}}\int_{\mathbb R}\frac{G_-(\omega')g_+(\omega')}{\omega'-{\rm i}}  {\rm d} \omega' =
\text{(perturbative terms)}
- \sum_{n\ge 1} \frac{{\rm e}^{-2B\xi_n} G_+({\rm i}\xi_n) {\rm i}\sigma_n^\pm g_n}{\xi_n-1} \\
-  h {\rm e}^{-2B} G_-({\rm i}) - \frac{m {\rm e}^{-B}}{2}\left(2BG_-({\rm i}) + {\rm i} G_-'({\rm i}\pm 0)\right),
\end{multline}
where
\begin{equation}
    g_n= -\frac{h}{\xi_n}+\frac{m{\rm e}^B}{2(\xi_n-1)}.
\end{equation}
For the particular value of $\kappa=1$, the non-integrated terms in \eqref{eps_eq} simplify to
\begin{equation}
\lim_{\kappa \rightarrow 1} \left[
\frac{h {\rm e}^{-2 B \kappa}}{\kappa G_+({\rm i} \kappa)}  -\frac{m{\rm e}^{B}}{2(\kappa-1)}\left(\frac{{\rm e}^{-2 B \kappa}}{G_+({\rm i}\kappa)}-\frac{{\rm e}^{-2 B}}{G_+({\rm i})}\right) \right] =  h \frac{{\rm e}^{-2 B}}{G_+({\rm i})} + \frac{m{\rm e}^{-B}}{2}\left(\frac{ 2B }{G_+({\rm i})}+\frac{ {\rm i} G_+'({\rm i})}{G^2_+({\rm i})}\right).
\end{equation}
Using the values $G_+({\rm i})$ and $G_-({\rm i})$ in \eqref{K+K+}, we can see that the above contributions change the sign in the real part of the two terms appearing in last line of \eqref{G-g+_computation}, while keeping their ambiguous imaginary part unmodified. More explicitly, when putting everything together, we find for the $\vartheta=\pi$ case
\begin{multline}
    \frac{\epsilon_+({\rm i})}{G_+({\rm i})} =
    \text{(perturbative terms)}
    - \sum_{n\ge 1} \frac{{\rm e}^{-2B\xi_n}{\rm i}\sigma_n^\pm v({\rm i}\xi_n)}{\xi_n-1}
    - \sum_{n\ge 1} \frac{{\rm e}^{-2B\xi_n} G_+({\rm i}\xi_n) {\rm i}\sigma_n^\pm g_n}{\xi_n-1}\\
    - \frac{{\rm e}^{-2B}  }{{\rm e}}v({\rm i})+h{\rm e}^{-2B}\sqrt{\frac{1}{2{\rm e}}}+\frac{m{\rm e}^{-B}}{2}\left( \sqrt{\frac{2}{{\rm e}}}   B-\frac{1  }{2\sqrt{2 {\rm e}}}\left(-1 +\gamma_E +\log 2\right)\right)\pm \frac{{\rm i}\pi m {\rm e}^{-B}}{4\sqrt{2{\rm e}}}.
    \label{eplus_expand}
\end{multline}
This contrasts with the $\vartheta=0$ case
\begin{multline}
    \frac{\epsilon_+({\rm i})}{G_+({\rm i})} =
    \text{(perturbative terms)}
    - \sum_{n\ge 1} \frac{{\rm e}^{-2B\xi_n}{\rm i}\sigma_n^\pm Q({\rm i}\xi_n)}{\xi_n-1}
    - \sum_{n\ge 1} \frac{{\rm e}^{-2B\xi_n} G_+({\rm i}\xi_n) {\rm i}\sigma_n^\pm g_n}{\xi_n-1}\\
    + \frac{{\rm e}^{-2B}  }{{\rm e}}Q({\rm i})-h{\rm e}^{-2B}\sqrt{\frac{1}{2{\rm e}}}-\frac{m{\rm e}^{-B}}{2}\left( \sqrt{\frac{2}{{\rm e}}}   B-\frac{1  }{2\sqrt{2 {\rm e}}}\left(-1 +\gamma_E +\log 2\right)\right)\pm \frac{{\rm i}\pi m {\rm e}^{-B}}{4\sqrt{2{\rm e}}},
    \label{eplus_theta0}
\end{multline}
where the perturbative terms are the same in both cases up to $Q\leftrightarrow v$.

Finally, in the $\kappa\rightarrow\infty$ limit, which we need to impose the boundary condition, the non-integrated terms with exponential dependency on $\kappa$ do not contribute.
From $\lim_{\kappa\rightarrow\infty} G_+({\rm i}\kappa)=1$, we find that the limit in the $\vartheta=\pi$ case is then given by
\begin{equation}
    \lim_{\kappa\rightarrow\infty} \kappa\epsilon_+({\rm i}\kappa) = \text{(perturbative terms)}+\sum_{n\ge 1} {\rm e}^{-2B\xi_n}{\rm i}\sigma_n^\pm \left(v({\rm i}\xi_n) + G_+({\rm i}\xi_n) g_n\right) - \frac{m{\rm e}^{-B}}{2\sqrt{2{\rm e}}},
    \label{bclim_pi}
\end{equation}
which contrasts with the $\vartheta=0$ case
\begin{equation}
    \lim_{\kappa\rightarrow\infty} \kappa\epsilon_+({\rm i}\kappa) = \text{(perturbative terms)}+\sum_{n\ge 1} {\rm e}^{-2B\xi_n}{\rm i}\sigma_n^\pm \left(Q({\rm i}\xi_n) + G_+({\rm i}\xi_n) g_n\right) + \frac{m{\rm e}^{-B}}{2\sqrt{2{\rm e}}}.
    \label{bclim_0}
\end{equation}

By comparing equations \eqref{veq_expand}, \eqref{eplus_expand} and \eqref{bclim_pi} with equations \eqref{Q_expand}, \eqref{eplus_theta0} and \eqref{bclim_0}, we see that
one can transform the integral equations for the problem with $\vartheta=\pi$ into the $\vartheta=0$ equations by swapping $Q\leftrightarrow v$ and ${\rm e}^{-2B}\rightarrow -{\rm e}^{-2B}$ in the
non-ambiguous terms, i.e. those not multiplied by ${\rm i} \sigma_n^\pm$. This is true
at the level of the equations and at leading order, but we can see that the solution
does not map as straightforwardly. For example, when solving for $v({\rm i}\xi)$ there is a term ${\rm e}^{-2B\xi_n}{\rm i} \sigma_n^\pm v({\rm i}\xi_n)$ which does not change sign under this map, but $v({\rm i}\xi_m)$ itself has a real non-perturbative correction of order ${\rm e}^{-2B}$ which does change sign. Furthermore, once we impose the boundary condition and re-express $B$ in terms of $h$, or alternatively $\tilde{\alpha}(h)$, the terms which swap and do not swap sign become further entangled. However, one can easily keep track of them by, for example, multiplying the appropriate terms by a factor $\cos(\vartheta)$ in the original equations.

From this prescription, we get the following result for the normalized energy density at $\vartheta=\pi$:
\begin{multline}
\frac{e}{\pi\rho^2} = \alpha+\frac{\alpha^2}{2} + \mathcal{O}\big(\alpha^3\big)
- \frac{32}{{\rm e}^2}\left( \frac{2}{\alpha } + \log(\alpha ) +3  - \gamma_E - 5\log(2) \pm \frac{ {\rm i}\pi}{2}
+\frac{\alpha}{2}
+\mathcal{O}\big(\alpha^2\big) \right)  {\rm e}^{-2/\alpha}\\
+ \frac{512}{{\rm e}^4}\left( \frac{1 \pm {\rm i}}{\alpha} + \mathcal{O}\big(\alpha^0\big) \right) {\rm e}^{-4/\alpha} + \mathcal{O}\big( {\rm e}^{-6/\alpha} \big),
\label{eq_normalized_energy_theta_pi}
\end{multline}
which is the expression we previously obtained for the $\vartheta=0$ case, in \eqref{eq_normalized_energy}, but changing the sign of the real part in the ${\rm e}^{-2/\alpha}$ correction. Let us recall that, in the presence of a
theta angle, an instanton configuration with topological charge ${\cal Q}$ in the $O(3)$ non-linear sigma model contributes
to a physical observable as
\begin{equation}
\exp\left\{ {\rm i} {\cal Q} \vartheta - {4 \pi |{\cal Q}| \over g^2(\mu)}\right\}.
\end{equation}
In particular, an instanton configuration with an odd value of ${\cal Q}$ will change sign as we go from $\vartheta=0$ to $\vartheta =\pi$. The change of sign of the real, exponentially small
correction of order ${\rm e}^{-2/\alpha}$ in (\ref{eq_normalized_energy_theta_pi})
is precisely what is expected from a configuration with ${\cal Q}=\pm 1$ (i.e. an (anti)-one-instanton). The non-ambiguous part of the correction of order ${\rm e}^{-4/\alpha}$, which does not change sign, would correspond to an instanton correction with ${\cal Q}=\pm 2$. More generally, the change of sign ${\rm e}^{-2B}\rightarrow -{\rm e}^{-2B}$ in the real contributions is the expected behavior of instanton corrections.

As was discussed in \cite{mmr-an}, one of the important roles of the formally imaginary ambiguous terms ${\rm i}\sigma_n^\pm$ is that they cancel the ambiguities originating in the resummation of the divergent perturbative part. It is thus important that they do not change sign between the two cases, since the perturbative part itself remains identical.  On the other hand, this also means that the terms that do swap sign are \textit{a priori} invisible to the perturbative series.

We want to remark that even if the residues ${\rm i}\sigma_n^\pm$ do not change sign, it does not follow that formally imaginary ambiguous terms do not change sign in the final result. As we mentioned earlier, terms like ${\rm e}^{-2B\xi_n}{\rm i} \sigma_n^\pm v({\rm i}\xi_n)$ will give rise to imaginary ambiguous terms that do change sign. This is still consistent with the cancellation of imaginary ambiguities in the Borel resummation of the whole trans-series. Indeed, if an exponential sector has a real part that changes sign, then the imaginary ambiguity arising from the Borel sum of that exponential sector has to change sign accordingly.

Let us note that our result for the change in the free energy as we vary the theta angle from
$\vartheta=0$ to $\vartheta=\pi$ suggests the following instanton-type dependence of the free energy on the theta angle $\vartheta$,
\begin{equation}
{\cal F}(h, \vartheta)=  -\frac{ h^2}{4 \pi}  \mathcal{F}_{(0)}(h) \mp {\rm i} \frac{m^2}{16}-\cos(\vartheta) \frac{ h^2}{4 \pi}  \mathcal{F}_{(1)}(h){\rm e}^{-2/\tilde\alpha} + \mathcal{O}\big({\rm e}^{-4/\tilde\alpha}\big),
\end{equation}
where $ \mathcal{F}_{(0)}(h)$ and $ \mathcal{F}_{(1)}(h)$ are given in (\ref{f1h}), and the term $\cos(\vartheta)$ comes from adding the contributions of the instantons with ${\cal Q}=\pm 1$. It is then natural to define an $h$-dependent topological susceptibility
\begin{equation}
\chi_t(h)= \left( {{\rm d}^2 {\cal F}(h, \vartheta) \over {\rm d} \vartheta^2} \right)_{\vartheta=0}.
\end{equation}
Therefore, we have the following asymptotic result as $\tilde \alpha \rightarrow 0$:
\begin{equation}
\label{chith}
\chi_t(h)\approx  \frac{ h^2}{4 \pi}  \left(\frac{64}{{\rm e}^2 \tilde\alpha^3} - \frac{32 (-\log(\tilde\alpha) -3 + \gamma_E +5\log (2))}{{\rm e}^2 \tilde\alpha^2}+O\big(\tilde\alpha^{-1}\big) \right){\rm e}^{-2/\tilde\alpha}.
\end{equation}
This is a genuine non-perturbative result for the $O(3)$ sigma model which may be testable in lattice computations (lattice calculations for the $O(3)$ sigma model with an $h$ field have been considered in e.g. \cite{bruckmann}).

\subsection{Testing the analytic results}
\label{sec_test_results_O3}

In the following section, we will give ample numerical evidence that the computations of sections \ref{sec_nonpert_theta_0} and \ref{sec_O3pi} led to the correct non-perturbative corrections for the $O(3)$ sigma model at $\vartheta=0$ and $\pi$. To this end, we compute numerically the normalized energy
density $e/(\pi\rho^2)$ and the associated $\alpha$ directly from the Bethe ansatz equation in their canonical formulation. The coupling $\alpha$ can be obtained from its definition in \eqref{eq_alpha_def}, together with the mass gap \eqref{mass-gap}, once $\rho$ has been computed.

For $\vartheta=0$, we solve numerically the integral equation \eqref{chi-ie} for $\chi(\theta)$. This can be done straightforwardly by discretizing the integral equation with a Gauss–Kronrod quadrature in the integration interval $[-B,B]$, and solving the resulting system of equations for $\chi(\theta)$. Then we
compute $\big(\alpha,e/(\pi\rho^2)\big)$ through the integrals of \eqref{erho_pi}, using the same quadrature rule. One can achieve very high precision in the computation by using discretizations of the integral equation with a relatively small number of points.

For $\vartheta=\pi$ the situation is very different. The integral equation \eqref{merTBA2} has an unbounded integration interval, and the kernel decays
slowly, like $1/x^2$. This leads to a substantial loss of precision. In order to achieve the accuracy required to detect exponentially small corrections, one needs a detailed analysis of the integral equation, and the development of appropriate numerical tools, which we explain in Appendix \ref{app-numerical-integral}.

Once we have obtained $\big(\alpha,e/(\pi\rho^2)\big)$ from the Bethe ansatz equations, we need to subtract the perturbative part,
so we can test our prediction for the first exponential correction. This requires some tools from the theory of resurgent asymptotics which we
now summarize briefly, see \cite{ss,mmlargen,mmbook,abs} for details and references.
Let us write the perturbative part of the normalized energy density as
\begin{equation}
\frac{e}{\pi\rho^2}\sim \alpha \varphi(\alpha), \qquad \varphi(\alpha) = \sum_{k\ge 0} e_k \alpha^k.
\end{equation}
We wrote the first two terms of this formal power series in \eqref{eq_normalized_energy} and \eqref{eq_normalized_energy_theta_pi}. In order to resum this series, we introduce its Borel transform
\begin{equation}
\widehat\varphi(\zeta) = \sum_{k\ge 0} { e_k  \over  k!} \zeta^k,
\end{equation}
which has singularities on the positive real axis, as shown in \cite{volin,mr-ren,mmr-an}. We also introduce
the (lateral) Borel resummations of $\varphi(\alpha)$ as
\begin{equation}
s_{\pm}(\varphi)(\alpha) = \frac{1}{\alpha} \int_{\mathcal{C}_\pm} \widehat\varphi(\zeta) {\rm e}^{-\zeta/\alpha} {\rm d} \zeta,
\end{equation}
where $\mathcal{C}_\pm$ is a straight line from 0 to $\infty$ slightly above or below the positive real line. The Borel sum
will have imaginary ambiguities due to the singularities in the Borel transform, and these cancel with the explicit imaginary
ambiguities appearing in the exponential corrections of \eqref{eq_normalized_energy} and
\eqref{eq_normalized_energy_theta_pi}. This cancellation was tested explicitly in \cite{mmr-an} for
the ambiguous term in the ${\rm e}^{-2/\alpha}$ correction and we checked with great precision
that this also happens for the ${\rm e}^{-4/\alpha}$ correction.

In \figref{fig_fullresultcompare}, we show the points $\big(\alpha,e/(\pi\rho^2)\big)$ that we obtained
by solving numerically the Bethe ansatz integral equations for the $O(3)$ sigma model, for both
$\vartheta=0$ and $\pi$. When $\vartheta=\pi$, we used the values $B=10/k$, with integers
$1 \le k \le 10$. When $\vartheta=0$, we looked for the values of $B$ that would result into
the same numerical $\alpha$ that we found for $\vartheta=\pi$ \footnote{Because the relation between $\alpha$ and $B$ includes exponential corrections, it is distinct between the two values of $\vartheta$.}. This can be done with a simple bisection method up to
desired precision. The dashed line corresponds to
$\alpha \Re\big[ s_{\pm}(\varphi)(\alpha) \big]$.  To obtain the Borel resummation of the perturbative series,
we used 70 coefficients, computed with the generalization of Volin's method \cite{volin, volin-thesis, mr-ren},
and we improved the numerical result with a conformal mapping, a strategy similar to the “Padé-Conformal-Borel”
method of \cite{costin-dunne-conformal}. The error estimates are extracted from the convergence of the Padé approximant, see \cite{mmr-an,dpmss}. One can see that $e/(\pi\rho^2)$ coalesces to the same
perturbative result for both values of the $\vartheta$ angle, when $\alpha$ is small (corresponding to large $B$).
However, as $\alpha$ grows (i.e. $B$ gets smaller), the two results start to differ, due to the different
exponential corrections to the perturbative expansion, see \figref{fig_fullresultcompare}.
\begin{figure}
\centering
\includegraphics[width=0.59\textwidth]{figures/compareO3pert.pdf}
\caption{Numerical evaluation of $\big(\alpha,e/(\pi\rho^2)\big)$ for the $O(3)$ sigma model with $\vartheta=0$ (black dots) and $\vartheta=\pi$ (red dots). The dashed line corresponds to the Borel resummation $\alpha \Re\big[ s_{\pm}(\varphi)(\alpha) \big]$ of the perturbative part, which is shared by both values of the $\vartheta$ angle.}
\label{fig_fullresultcompare}
\end{figure}

We can test the ${\rm e}^{-2/\alpha}$ correction in our results by considering the difference between
$e/(\pi\rho^2)$ and the real part of the Borel sum. According to \eqref{eq_normalized_energy}
and \eqref{eq_normalized_energy_theta_pi}, this difference should have the asymptotic expansion
\begin{multline}
\frac{e}{\pi\rho^2} - \alpha\Re\big[ s_{\pm}(\varphi)(\alpha) \big] \sim \pm \frac{32}{{\rm e}^2}\biggl(\frac{2}{\alpha } + \log(\alpha ) +3  - \gamma_E - 5\log(2)\\
+ \frac{\alpha}{2} + \mathcal{O}\big(\alpha^2\big) \biggr)  {\rm e}^{-2/\alpha}
+ \mathcal{O}\big( {\rm e}^{-4/\alpha} \big),
\label{erho-exp-correction}
\end{multline}
where the $+$ and $-$ sign correspond to $\vartheta=0$ and $\pi$, respectively. We have also included the term of order $\alpha$ that was calculated in \cite{bbh}. The asymptotic expansion for $\vartheta=0$  was already tested in \cite{mmr-an,bbh}, so we focus our attention on the $\vartheta=\pi$ case.

In \figref{fig_e_rho_exp_O3}, we compare the numerical value of
$e/(\pi\rho^2) - \alpha\Re\big[ s_{\pm}(\varphi)(\alpha) \big]$ against the asymptotic
expansion of \eqref{erho-exp-correction}, for $\vartheta=\pi$. The same information is shown
in Table \ref{tab_e_rho_exp_O3}. In general, we find a very good agreement. For large values of $\alpha$ there is a slight
disagreement that can be accounted for missing terms in the asymptotic expansion. For the points with small $\alpha$, there is also a
slight disagreement (in absolute terms) that can only be appreciated in the plot on the right in
\figref{fig_e_rho_exp_O3}\footnote{The point with lowest $\alpha$ is completely off in this plot, and it does not show inside the displayed range, although its estimated error is similarly very large. Small numerical errors become large when multiplying by ${\rm e}^{2/\alpha}$, with
$\alpha$ small, thus very high numerical precision is required.}. In this case, the error originates in the numerical computation of $e$ and $\rho$,
and it could be reduced by allocating more computer time into the integral equations of the Bethe ansatz. The error bars
have been estimated from the convergence of $\rho$ when increasing the order of the integral equation quadrature, as explained at the end of Appendix \ref{app-numerical-integral}.
\begin{figure}
\centering
\includegraphics[width=0.43\textwidth]{figures/e_rho_exp_O3.pdf}
\hspace{0.065\textwidth}
\includegraphics[width=0.43\textwidth]{figures/e_rho_exp_norm_O3.pdf}
\caption{In the figure on the left we plot the difference between the normalized energy density $e/(\pi\rho^2)$ and the real
part of the Borel resummation of the perturbative series, for the $O(3)$ sigma model at $\vartheta=\pi$. The dots are the numerical calculations of the difference. The two lines are the theoretical prediction, given by the asymptotic expansion in the r.h.s. of \eqref{erho-exp-correction}, including terms up to order $\alpha^0$ (dashed) and $\alpha$ (solid). In the figure on the right, the result has been normalized by dividing by ${\rm e}^{-2/\alpha}\, 64/({\rm e}^2\alpha)$.}
\label{fig_e_rho_exp_O3}
\end{figure}

\begin{table}
\begin{equation*}
\arraycolsep=2.5mm
\begin{array}{lll}
\toprule
 \alpha & \text{Numerical result} & \text{Asymptotics}\\
\midrule
 0.082486 & (0.5\! \pm \! 2.6)\times 10^{-6} & -2.65089\times 10^{-9}\\
 0.146461 & -0.00005(4\! \pm \! 5) & -0.000055\\
 0.200073 & -0.00146(8\! \pm \! 7) & -0.001469\\
 0.246564 & -0.00752(2\! \pm \! 9) & -0.007526\\
 0.287605 & -0.0198(57\! \pm \! 10) & -0.019884\\
 0.324195 & -0.0376(11\! \pm \! 11) & -0.037725\\
 0.357020 & -0.0589(40\! \pm \! 12) & -0.059259\\
 0.386608 & -0.0820(84\! \pm \! 13) & -0.082778\\
 0.413390 & -0.1057(27\! \pm \! 14) & -0.106994\\
 0.437725 & -0.1290(03\! \pm \! 15) & -0.131052\\
\bottomrule
\end{array}
\end{equation*}
\caption{Difference between the normalized energy density $e/(\pi\rho^2)$ and the real
part of the Borel resummation of the perturbative series, for the $O(3)$ sigma model at $\vartheta=\pi$. We compare the numerical value obtained from the Bethe ansatz (2nd column), to the asymptotic expansion in the r.h.s. of \eqref{erho-exp-correction} (3rd column).}
\label{tab_e_rho_exp_O3}
\end{table}

We also tested the ${\rm e}^{-4/\alpha}$ correction in \eqref{eq_normalized_energy} and \eqref{eq_normalized_energy_theta_pi}.
As we do not have enough terms in the $\alpha$ expansion accompanying ${\rm e}^{-2/\alpha}$, it
is not possible to subtract the whole exponential correction and still be left with enough precision
to test the subleading exponential correction. However, one way to circumvent this issue is to consider instead
the sum of $e/(\pi\rho^2)$ at $\vartheta=0$ and $\vartheta=\pi$. According to \eqref{eq_normalized_energy}
and \eqref{eq_normalized_energy_theta_pi}, the correction of order ${\rm e}^{-2/\alpha}$ cancels in the sum, leading to
\begin{equation}
\frac{1}{2\pi}\left[ \frac{e}{\rho^2}(\alpha,0) + \frac{e}{\rho^2}(\alpha,\pi) \right] -  \alpha\Re\big[ s_{\pm}(\varphi)(\alpha) \big] \sim \left( \frac{512}{{\rm e}^4 \alpha} + \mathcal{O}\big(\alpha^0\big) \right) {\rm e}^{-4/\alpha} + \mathcal{O}\big( {\rm e}^{-8/\alpha} \big).
\label{erho-exp2-correction}
\end{equation}
We test this expression in \figref{fig_e_rho_exp2_O3}, finding again good agreement.
\begin{figure}
\centering
\includegraphics[width=0.59\textwidth]{figures/secondsectorO3.pdf}
\caption{Difference between $\big[ e/\rho^2(\alpha,0) + e/\rho^2(\alpha,\pi) \big]/(2\pi)$ and the real part of the Borel resummation of the perturbative series, for the $O(3)$ sigma model. The dots are the numerical calculations of the difference. The dashed line is the theoretical prediction, given by the asymptotic expansion in the r.h.s. of \eqref{erho-exp2-correction}.}
\label{fig_e_rho_exp2_O3}
\end{figure}

\sectiono{Fendley's integrable models}
\label{fendley-sec}

\subsection{General aspects}

 Due to the integrability of the $O(3)$ sigma model at $\vartheta =0$ and $\vartheta= \pi$, we have been able to study the effects
 of the theta angle on the non-perturbative effects appearing in the free energy. It is natural to ask whether this feature can be
 generalized to other models. This was answered in the affirmative in \cite{fendley}. The coset sigma models
$SU(N)/SO(N)$ and $O(2P)/O(P) \times O(P)$ have
$\pi_1(G/H)= {\mathbb Z}_2$ (for $N, P>2$), so there is a ${\mathbb Z}_2$ topological charge and one can introduce a
discrete theta angle which can only take the values $\vartheta=0, \pi$. Fendley conjectured that the models
are integrable and provided conjectural expressions for the $S$ matrices, for both values of $\vartheta$ (see also \cite{fendley2}).

Let us now review some useful results of \cite{fendley}. We start by reminding some basic facts about coset models
based on symmetric spaces $G/H$. Their action can be written as
\begin{equation}
S={1\over g_0^2} \int {\rm d}^2 x \Tr \left( \partial_\mu U^{-1} (x) \partial_\mu U(x) \right),
\end{equation}
where $U$ is a matrix field taking values in the Lie group $G$ and satisfying some additional constraints.
In the case of the $SU(N)/SO(N)$ coset model,
$U$ is a symmetric $SU(N)$ matrix. In the case of the $O(2P)/O(P) \times O(P)$ coset model, $U$ is a symmetric and traceless $O(2P)$ matrix.
In the presence of a chemical potential $h$ coupled to a conserved charge $Q$, the action is modified as \cite{pcf}
\begin{equation}
S={1\over g_0^2} \int {\rm d}^2 x \Tr \left(\overline D_\mu U^{-1} (x) D_\mu U(x) \right),
\end{equation}
where
\begin{equation}
\begin{aligned}
D_\mu U &= \partial_\mu U - h \delta_{\mu 0} (Q U + U Q), \\
\overline D_\mu U^{-1}&=  \partial_\mu U^{-1} + h \delta_{\mu 0} (U^{-1} Q + Q U^{-1}),
\end{aligned}
\end{equation}
and the charge $Q$ is an element in the Cartan subalgebra $\mathfrak{h}$ of the Lie algebra $\mathfrak{g}$ of $G$.
It is convenient to represent $Q$ by its dual, which is
a vector $\boldsymbol{q}\in \Lambda_{\rm w} \otimes {\mathbb R}$. Here, $\Lambda_{\rm w}$ is the weight lattice. One can then calculate
the free energy ${\cal F}(h)$ at one loop and obtain \cite{pcf,hollo-lie}
\begin{equation}
\label{fh-lie}
{\cal F}(h)= -{4 h^2 \over g^2(h)}  \boldsymbol{q}^2-{h^2 \over 2 \pi} \sum_{{\boldsymbol{\alpha}} \in \Delta_+} (\boldsymbol{q} \cdot {\boldsymbol{\alpha}})^2 \left\{
\log \left( |\boldsymbol{q}\cdot{\boldsymbol{\alpha}} | \right)-{1\over 2} \right\}+\cdots,
\end{equation}
where $\Delta_+$ is the set of positive roots of $\mathfrak{g}$ and $g^2(h)$ is the running coupling constant at the scale set by $h$.

In \cite{fendley}, Fendley considered the free energy ${\cal F}(h)$ for the two coset models $SU(N)/SO(N)$
and $O(2P)/O(P) \times O(P)$, with particular choices of the charge $\boldsymbol{q}$. By comparing the perturbative result (\ref{fh-lie}) to the Bethe ansatz
calculation of ${\cal F}(h)$, he checked his proposals for the $S$ matrices of the theory. Let us first consider the results for the $SU(N)/SO(N)$ coset sigma model. This model is asymptotically free, and in the convention (\ref{betaf}) for the beta function, we have
\begin{equation}
\label{beta-coset}
\beta_0= {1\over 8 \pi \Delta}, \qquad \beta_1= {1\over 128 \pi^2} {1+2 \Delta \over  \Delta^2},
\end{equation}
where
\begin{equation}
\Delta={1\over N}.
\end{equation}
We now consider the coupling to a conserved charge. When $G=SU(N)$, the charge $\boldsymbol{q}$ can be written as a vector in ${\mathbb R}^N$,
\begin{equation}
\boldsymbol{q}=(q_1, \cdots, q_N), \qquad \sum_{j=1}^n q_j=0,
\end{equation}
and the free energy (\ref{fh-lie}) is then given by
\begin{equation}
\label{fsun}
{\cal F}(h)=-{4 h^2 \over g^2(h)}  \sum_{j=1}^N q_j^2 -{h^2 \over 2 \pi} \sum_{i<j} \left(q_i-q_j\right)^2 \left\{ \log\left( |q_i - q_j|  \right)-{1\over 2} \right\}+\cdots.
\end{equation}
In \cite{fendley}, the following choice is made for the charge,
\begin{equation}
\label{qcharge}
\boldsymbol{q}= \mathsf{r}(1, -1, 0, \cdots,0),
\end{equation}
where $\mathsf{r}$ is a normalization factor which is chosen to match the Bethe ansatz calculation. Let us note that with this choice of charge there
are two different types of particles in the ground state. By using the exact $S$-matrices proposed in \cite{fendley}, one finds that the
free energy can be calculated from the integral equation (\ref{eps-ie}) with the kernel
\begin{equation}
\widetilde K(\omega) = {\rm e}^{-\pi\Delta|\omega|} \frac{2\cosh\big( \tfrac{1}{2}(1-2\Delta)\pi\omega \big)
\sinh(2\pi\Delta|\omega|)}{\cosh\big(\tfrac{1}{2}\pi\omega\big)}.
\label{eq_kernel}
\end{equation}
The Wiener--Hopf factorization of this kernel is given by\footnote{This corrects some minor errors in the corresponding equation in \cite{fendley}.}
\begin{multline}
G_+(\omega) =\frac{ \exp \left({\rm i} \omega  \left(-\Delta +2 \Delta  \log (2 \Delta )+\left(\frac{1}{2}-\Delta \right) \log \left(\frac{1}{2}-\Delta \right)+\frac{\log (2)}{2}\right)+{\rm i} \Delta  \omega  \log (-{\rm i} \omega )\right)}{2 \sqrt{-\pi {\rm i} \Delta\omega }}\\
\times\frac{\Gamma (1-2 {\rm i} \Delta  \omega ) \Gamma \left(\frac{1}{2}-{\rm i} \left(\frac{1}{2}-\Delta \right) \omega \right) }{ \Gamma \left(\frac{1}{2}-\frac{{\rm i} \omega }{2}\right)}.
\label{Gplus}
\end{multline}
The perturbative free energy can be computed from the Bethe ansatz at one loop by using (\ref{ftba}). Since there are two particles in the ground state,  we have
to include an additional factor of $2$ in (\ref{ftba}), and by comparing to (\ref{fsun}) we find that
the normalization of the charge in (\ref{qcharge}) is $\mathsf{r}^2=1/8$
(a similar normalization was found in \cite{babichenko}). Agreement of the two expressions requires that \cite{fnw1}
\begin{equation}
\label{beta-cond}
{\beta_1 \over \beta_0^2}=a+{1\over 2},
\end{equation}
which holds in this model since $a= \Delta$,
and tests the correctness of the $S$-matrix ansatz proposed in \cite{fendley}.
From these results one can obtain in addition the mass gap, which is given by
\begin{equation}
{m \over\Lambda}= {\sqrt { \pi}}  {2^{3\Delta+2 } {\rm e}^{-{1\over 2}-\Delta } \over \Gamma(1- \Delta) \Gamma(1+ 2 \Delta)}.
\end{equation}

The $O(2P)/O(P) \times O(P)$ coset sigma model has the same values for
the coefficients of the beta function given in (\ref{beta-coset}), with the only difference
that now
\begin{equation}
\label{delta-O}
\Delta={1\over 2(P-1)}.
\end{equation}
In \cite{fendley} the following charge is considered
\begin{equation}
\boldsymbol{q}= \mathsf{r} \boldsymbol{\lambda}_1,
\end{equation}
where $\boldsymbol{\lambda}_1$ is the highest weight of the fundamental representation of
$O(2P)$, and $\mathsf{r}$ is a normalization constant. The ground state is then populated by a single species of particles, see \cite{fendley} for details.
For this choice of the charge, the kernel appearing in the Bethe ansatz
has the same form than (\ref{eq_kernel}), where $\Delta$ is now given in (\ref{delta-O}). This makes it possible to analyze both models
simultaneously. By matching the perturbative and the Bethe ansatz answers for the free energy
one checks (\ref{beta-cond}). We find $\mathsf{r}^2=1/8$, as well as the following value for the mass gap,
\begin{equation}
{m \over\Lambda}= {\sqrt { \pi}}  {2^{5 \Delta+2 } {\rm e}^{-{1\over 2}-\Delta } \over \Gamma(1- \Delta) \Gamma(1+ 2 \Delta)}.
\end{equation}

\subsection{The models at \texorpdfstring{$\vartheta=0$}{theta = 0}}
\label{sec_fen0}

At $\vartheta=0$, the integral equations for these models are structurally identical to the ``bosonic models''  studied in \cite{mr-ren,mmr-an} (which include the $O(N)$ sigma model, see also \cite{bbhv22}), and so we can easily apply the perturbative and non-perturbative analysis therein. From this section on,
we will consider the observables as defined per excited particle species. This changes nothing in the case of $O(2P)/O(P)\times O(P)$, but in the case of
$SU(N)/SO(N)$ there are two particle types in the ground state that obey the same Bethe ansatz equation and so the total free energy is twice the free energy per particle, which we present, and so forth.
To express these results, we introduce a coupling $\alpha$, analogous to \eqref{eq_alpha_def},
\begin{equation}
\frac{1}{\alpha }+(\xi -1) \log (\alpha )=\log \left(\frac{C \rho }{m}\right),
\label{alpha_def}
\end{equation}
where
\begin{equation}
    \xi = \frac{\beta_1}{\beta_0^2},\quad C = \frac{m}{\Lambda}\frac{1}{\mathfrak{c}\beta_0},
\end{equation}
using the convention \eqref{betaf} (which is distinct from the convention used in \cite{mr-ren,mmr-an}).
In both families of coset models, we choose the arbitrary constant $\mathfrak{c}$ such that
\begin{equation}
C=\frac{\pi ^{3/2} 2^{5 \Delta +\frac{7}{2}} {\rm e}^{-\Delta -1/2} \Delta }{\Gamma (1-\Delta ) \Gamma (2 \Delta +1)},
\end{equation}
since with this choice the perturbative coefficients take the simplest possible form.

Using \eqref{Gplus}, we can apply the generalised Volin's method from \cite{mr-ren,volin} to obtain the perturbative expansion of the ground state energy to high order. Following the procedure for ``bosonic models'', we find
\begin{multline}
\frac{e}{\rho^2} = \frac{\alpha}{4\pi\Delta}\bigg[
1+\frac{\alpha }{2}+\alpha ^2 \left(\frac{1}{4}+\frac{\Delta }{2}\right)
+\alpha ^3 \bigg(\frac{5}{16}-\frac{3}{32} (7 \zeta (3)-10) \Delta +\frac{1}{16} (21 \zeta (3)+8) \Delta ^2
\\
+\frac{1}{4} (-3 \zeta (3)-1) \Delta ^3\bigg)
+\alpha ^4 \bigg(\frac{53}{96}+\left(\frac{197}{96}-\frac{63 \zeta (3)}{64}\right) \Delta +\left(\frac{97}{48}-\frac{35 \zeta (3)}{16}\right) \Delta ^2
\\
+\left(\frac{115 \zeta (3)}{16}+\frac{1}{8}\right) \Delta ^3+\left(-\frac{19 \zeta (3)}{4}-\frac{1}{4}\right) \Delta ^4\bigg)
+{\cal O}\left(\alpha ^5\right)
\bigg].
\label{epert0}
\end{multline}
An interesting property of this series is that the leading part of the coefficients in the large $N$ limit, i.e.
$\Delta\rightarrow 0$, is identical to the leading part of the same coefficients in the large $N$ limit of the
PCF model, see \cite{mmr-an,dpmss}. We also note that the case $\Delta=1/2$ is special and it reduces to the
$O(3)$ sigma model. However, the subsequent formulae do not apply to that case.

Let us now consider the non-perturbative corrections. According to the general results of \cite{mmr-an},
the trans-series for the normalized ground state energy is given by the formula
\begin{equation}
\label{egs_trans}
\frac{e}{\rho^2}=\frac{\alpha}{k^2}\left\{\varphi(\alpha)+{\cal C}_0^\pm{\rm e}^{-\frac{2}{\alpha}}\alpha^{1-2\xi}+{\cal C}_1^\pm\left({\rm e}^{-\frac{2}{\alpha}}\alpha^{1-2\xi}\right)^{\xi_1}+\cdots \right\}
\end{equation}
where $\xi_1$ is the position of the first pole of $G_-({\rm i}\xi)$. Using the definition \eqref{alpha_def}, the constants are
\begin{align}
\label{C0_generic}
{\cal C}_0^\pm &= -\frac{C^2 G_+({\rm i}) k^2 }{4 \pi}{\rm i} G'_-({\rm i}\pm 0),\\
{\cal C}_1^\pm &=\frac{ 2{\rm i} \sigma^\pm _1 }{\left(\xi _1-1\right)^2 \xi _1}\left(\frac{C^2 G^2_+({\rm i}) k^2}{8\pi}\right)^{\xi _1},
\label{C1generic}
\end{align}
where $\sigma_1^\pm$ is the residue of $\sigma(\omega)$ at the pole $\xi_1$.
The free energy can be similarly written as
\begin{equation}
\label{fh-bos}
\begin{aligned}
{\cal F} (h) &= - \frac{k^2 h^2}{4\tilde\alpha}\Biggl\{ 1
-
{\rm i} \sigma^\pm_1\left({\rm e}^{-\frac{2}{\tilde{\alpha}}}\tilde{\alpha}^{1-2 \xi }\right)^{\xi _1}\!
\frac{2}{\left(\xi _1-1\right)^2 \xi _1}\left(\frac{G^2_+({\rm i})}{2 \pi  k^2}
\left(\frac{m}{\Lambda}\right)^{2}\right)^{\xi_1}
+ \cdots \Biggr\}\\
&\qquad - \frac{m^2}{4\pi} {\rm i} G_+({\rm i}) G'_-({\rm i}\pm 0),
\end{aligned}
\end{equation}
where $\tilde\alpha$ is given by
\begin{equation}
\frac{1}{\tilde{\alpha}}+\xi \log\tilde\alpha = \log\left( \frac{h}{\Lambda} \right).
\end{equation}
We can now specialize these general formulae to the Fendley's models. In these models, the poles of $\sigma(\omega)$ for any $\Delta$ are located
at
\begin{equation}
\frac{\ell}{1-2\Delta},\quad \frac{\ell'}{2\Delta},\quad \ell,\ell' \in {\mathbb N}.
\label{Gplus_poles}
\end{equation}
The residues of $\sigma(\omega)$ associated to the first family formally contribute as terms with both real and (ambiguous)
imaginary parts, while the second are purely real.
For $\Delta < 1/4$ the leading pole occurs at
\begin{equation}
\xi_1 = \frac{1}{1-2\Delta},
\label{xi1_fen}
\end{equation}
and consequently the coefficients (\ref{C0_generic}), (\ref{C1generic}) read
\begin{equation}
\begin{aligned}
{\cal C}_0^\pm(0) &=
\bigl[\sin (\pi  \Delta )\mp{\rm i}  \cos (\pi  \Delta )\bigr]
\left(\frac{64}{{\rm e}^2}\right)^{\Delta }\frac{ \Gamma \left(\frac{1}{2}-\Delta \right) }{\left({\rm e} \Delta \right) \Gamma \left(\Delta +\frac{1}{2}\right)},\\
{\cal C}_1^\pm(0) &=
\left[-\sin \biggl(\frac{\pi  \Delta }{1-2 \Delta }\biggr)\pm{\rm i} \cos \biggl(\frac{\pi  \Delta }{1-2 \Delta }\biggr)\right]
 \left(\frac{1024^{\Delta }}{4 {\rm e}^2 (1-2 \Delta )}\right)^{\frac{1}{1-2 \Delta }}
\frac{2 {\rm e} \left(1-2\Delta\right) \Gamma \bigl(\frac{1-4 \Delta }{2-4 \Delta }\bigr) }{ \Delta \Gamma \bigl(\frac{3-4 \Delta }{2-4 \Delta }\bigr)}.
\end{aligned}
\label{fenCpm_0}
\end{equation}
The argument in the l.h.s. is to remind that these are the coefficients at $\vartheta=0$. By using the asymptotic behavior of \eqref{epert0} and the numerical solution of the integral equation \eqref{chi-ie},
we have tested the real and imaginary parts of ${\cal C}_{0,1}^\pm$ with great precision for several values of $\Delta$, as in \cite{mmr-an}. The cases of $\Delta=1/3$, $1/4$ are conceptually similar but have some technical peculiarities, although formulae \eqref{C0_generic} and \eqref{C1generic} still apply.

Let us finally note that, following \cite{zamo-mass,foz,mmr-an}, the real part of the $h$-independent term in \eqref{fh-bos}
should give the free energy at $h=0$, and we find
\begin{equation}
F(0) = - \frac{m^2}{16\sin(2\pi\Delta)}.
\end{equation}

\subsection{The models at \texorpdfstring{$\vartheta =\pi$}{theta = pi}}

Since these models remain integrable in the presence of the $\vartheta=\pi$ angle, one can analyse their Bethe ansatze, similar to what we did for the $O(3)$ sigma model. As we will show, the resulting picture is also the same: the perturbative part does not change, non-perturbative formally imaginary ambiguous terms are also unchanged, but the non-perturbative formally real terms change sign. However, this change
is organized differently. In the $O(3)$ sigma model, the residues ${\rm i} \sigma_n^\pm$ were purely imaginary and the real non-perturbative terms came from explicit driving terms. In Fendley's sigma models, the latter do not appear, but the residues have both a real and imaginary part. We will thus first show that, when compared to the $\vartheta=0$ case, only the residues change, and then we will show that this change is only in the sign of the real part.

Up to equation \eqref{completeeq_genform}, the discussion in section \ref{sec_O3pi} also applies to Fendley's coset models. In this case, the kernels are given by \cite{fendley}
\begin{equation}
\phi_1(\omega) = \frac{1-{\rm e}^{2\pi\Delta|\omega|}}{1+{\rm e}^{\pi|\omega|}},
\qquad
\phi_2(\omega) = \frac{1+{\rm e}^{\pi(1-2 \Delta ) |\omega| }}{1+{\rm e}^{\pi|\omega|}}.
\end{equation}
and we easily verify \eqref{phitoK}, where now $\tilde K(\omega)$ is the kernel in \eqref{eq_kernel} for the Fendley's coset models at $\vartheta=0$. However, we no longer have the simplification of \eqref{phi2phi1_simplification} that we found for $O(3)$. Instead, we find
\begin{equation}
\frac{\phi_2(\omega)}{1-\phi_1(\omega)}= {\rm e}^{-2\pi\Delta|\omega|}.
\label{phi2phi1_fendley}
\end{equation}
In order to extend this function to the complex plane, we define
\begin{equation}
\frac{\phi_2(\omega)}{1-\phi_1(\omega)}= 1-\frac{1}{d_+(\omega)d_-(\omega)},
\label{eq_ddef}
\end{equation}
where
\begin{equation}
d_+(\omega) = \frac{{\rm e}^{-{\rm i}\Delta\omega[1-\log(-{\rm i}\Delta\omega)]}}{\sqrt{-2\pi{\rm i}\Delta\omega}} \Gamma(1-{\rm i}\Delta\omega),
\label{eq_dformula}
\end{equation}
and $d_-(\omega)=d_+(-\omega)$. An important ingredient for the analysis of the non-perturbative corrections at $\vartheta=\pi$ is the following result:
\begin{equation}
\frac{G_-({\rm i}\xi)}{d_+({\rm i}\xi)d_-({\rm i}\xi)} = 2 \Re\big(G_-({\rm i}\xi)\big), \qquad \xi>0.
\label{eq_d+d-_mod}
\end{equation}
This relation is easy to prove by observing that
\begin{equation}
d_+({\rm i} \xi \mp  0) d_-({\rm i} \xi \mp  0) = \frac{{\rm e}^{\pm{\rm i}\pi\big(\Delta\xi-\tfrac{1}{2}\big)}}{2\sin\big(\pi \Delta\xi \big)}, \qquad \xi>0.
\end{equation}
At the same time, we can write \eqref{Gplus} as
\begin{equation}
G_-({\rm i}\xi\mp 0) = r(\xi) {\rm e}^{\pm{\rm i}\pi\big(\Delta\xi-\tfrac{1}{2}\big)},
\end{equation}
where $r(\xi)$ is a real function. Then it follows that
\begin{equation}
\frac{G_-({\rm i}\xi)}{d_+({\rm i}\xi)d_-({\rm i}\xi)} = 2\sin\big(\pi\Delta\xi\big) r(\xi), \qquad \xi>0,
\end{equation}
from where we conclude the result of \eqref{eq_d+d-_mod}.

Following a similar script to section \ref{sec_O3pi}, we take the $(+)$-projection of \eqref{completeeq_genform},
\begin{multline}
v= -\left[{\rm e}^{-2{\rm i} B \omega}\left\{\sigma\frac{\phi_2}{1-\phi_1}v\right\}(-\omega)\right]_+ \\
+ \left[\frac{{\rm i} h G_+}{\omega-{\rm i} 0}\left(1-{\rm e}^{-2{\rm i} B \omega} \frac{\phi_2}{1-\phi_1}\right)\right]_+
- \frac{{\rm i} m {\rm e}^B}{2} \left[G_+\left(\frac{1}{\omega-{\rm i}}-\frac{\phi_2}{1-\phi_1}\frac{{\rm e}^{-2{\rm i} B \omega}}{\omega+{\rm i}}\right)\right]_+.
\label{vnotequal}
\end{multline}
We focus first on the integrals with $v$ in \eqref{vnotequal}. The function
\begin{equation}
\sigma(\omega) \frac{\phi_2(\omega)}{1-\phi_1(\omega)} = \sigma(\omega) - \frac{G_-(\omega)}{G_+(\omega)d_+(\omega)d_-(\omega)}
\label{atilde_Fendley}
\end{equation}
is the analogous of \eqref{atilde} in the manipulations of the $O(3)$ case. Now, the extra term in \eqref{atilde_Fendley} does have singularities in the upper half plane. However, it still has no branch cut, given that the ambiguities cancel in \eqref{eq_d+d-_mod}. In addition, the singularities of the extra term are located at ${\rm i}\xi_n$, the same positions as $\sigma(\omega)$, with residues $\tilde\sigma_n$. Since this function is not discontinuous, these residues are formally unambiguous. We can then deform the contour around the positive imaginary axis as usual, picking up the discontinuity of $\sigma(\omega)$ and the residues of both terms in \eqref{atilde_Fendley}:
\begin{equation}
\left[{\rm e}^{-2{\rm i} B \omega}\left\{\sigma \frac{\phi_2}{1-\phi_1}v\right\}(-\omega)\right]_+ =
\frac{1}{2\pi{\rm i}}\int_{{\cal C}_\pm} \frac{{\rm e}^{-2  B\xi'}\delta\sigma({\rm i}\xi')v({\rm i}\xi')}{\xi'-{\rm i}\omega} {\rm d} \xi' - \sum_{n\ge 1} \frac{{\rm e}^{-2B\xi_n}({\rm i}\sigma^\pm_n-{\rm i}\tilde\sigma_n) v({\rm i}\xi_n)}{\xi_n-{\rm i}\omega}.
\end{equation}
As for the driving terms, similar manipulations give the terms proportional to $m$ as
\begin{multline}
\left[G_+\left(g^m_-+{\rm e}^{-2{\rm i} B \omega}\frac{\phi_2}{1-\phi_1}g_+^m
\right)\right]_+
=
- \frac{{\rm i} m {\rm e}^B}{2} \frac{1}{2\pi{\rm i}} \int_{{\cal C}_\pm}\frac{{\rm e}^{-2B\xi'}\delta G_-({\rm i}\xi')}{({\rm i}\xi'+\omega)(\xi'-1)} {\rm d}\xi'
\\+\frac{{\rm i} m {\rm e}^B}{2} \frac{G_+(\omega)-G_+({\rm i})}{\omega-{\rm i}}
 + \frac{{\rm i} m {\rm e}^B}{2}\sum_{n\ge 1} \frac{{\rm e}^{-2B\xi_n }G_+({\rm i}\xi_n)({\rm i}\sigma^\pm_n-{\rm i}\tilde\sigma_n)}{({\rm i}\xi_n+\omega)(\xi_n-1)},
\end{multline}
and, with some additional work\footnote{The integral with the discontinuity is naively divergent due to the $1/\omega$ factor, however this can be fixed with the mouse-hole contour prescription described in Appendix A of \cite{mmr-an}.}, the terms proportional to $h$ can also be written as
\begin{multline}
\left[\frac{{\rm i} h G_+}{\omega-{\rm i} 0}\left(1-{\rm e}^{-2{\rm i} B \omega}\frac{\phi_2}{1-\phi_1}\right)\right]_+
=  {\rm i} h \frac{G_+(\omega)}{\omega} - \frac{{\rm i} h}{2\pi{\rm i}}\int_{{\cal C}_\pm} \frac{{\rm e}^{-2B\xi'}\delta G_-({\rm i}\xi')}{\xi'({\rm i}\xi'+\omega)} {\rm d}\xi'
\\
+ {\rm i} h \sum_{n\ge 1} \frac{{\rm e}^{-2B\xi_n}G_+({\rm i}\xi_n)({\rm i}\sigma^\pm_n-{\rm i}\tilde\sigma_n)}{\xi_n({\rm i}\xi_n+\omega)}.
\end{multline}
The $\vartheta=0$ case is recovered by removing the terms with ${\rm i}\tilde\sigma_n$. Thus, we can transform the $\vartheta=0$ equation for $Q$ into the $\vartheta=\pi$ equation for $v$ simply by swapping $\{Q,{\rm i}\sigma_n^\pm\}\leftrightarrow\{v,{\rm i}\sigma_n^\pm-{\rm i}\tilde\sigma_n\}$.

The equation for $\epsilon_-$, which follows from $(-)$-projection of \eqref{completeeq_genform}, is
\begin{multline}
\frac{\epsilon_-(\omega)}{G_-(\omega)} = \left[{\rm e}^{-2{\rm i} B \omega}\left\{\sigma v\right\}(-\omega)\right]_- + \left[G_+\left(g^m_-+{\rm e}^{-2{\rm i} B \omega}g^m_+\right)\right]_-
-{\rm i} h \left[G_+\frac{1-{\rm e}^{-2{\rm i} B \omega}}{\omega-{\rm i} 0}\right]_-\\
-\left[{\rm e}^{-2{\rm i} B \omega}\left\{\frac{G_- v}{G_+d_+d_-}\right\}(-\omega)\right]_-
-\left[{\rm e}^{-2{\rm i} B \omega}\left\{\frac{G_-(G_+g_-^m)}{G_+d_+d_-}\right\}(-\omega)\right]_-
-\left[{\rm e}^{-2{\rm i} B \omega}\frac{G_+}{G_-d_+d_-}\frac{{\rm i} h G_-}{\omega-{\rm i} 0}\right]_-.
\label{epsilon_minus}
\end{multline}
The terms in the first line are identical to the $\vartheta=0$ case, while those in the second are particular to $\vartheta=\pi$. We can write them in the generic form
\begin{equation}
\left[{\rm e}^{-2{\rm i} B \omega}\left\{\frac{G_- f}{G_+d_+d_-}\right\}(-\omega)\right]_- = \frac{1}{2\pi{\rm i}}\int_{\mathbb R} \frac{{\rm e}^{2{\rm i} B\omega'}G_-(\omega')}{G_+(\omega')d_+(\omega')d_-(\omega')}\frac{f(\omega)}{\omega'+\omega+{\rm i} 0} {\rm d} \omega',
\end{equation}
where
\begin{equation}
    f(\omega) = v(\omega) - {\rm i} h \frac{G_+(\omega)}{\omega+{\rm i} 0} + \frac{{\rm i} m{\rm e}^B}{2}\frac{G_+(\omega)}{\omega-{\rm i}}
\end{equation}
is a function analytic in the upper half plane except at $\omega={\rm i}$.  Since $G_-({\rm i})=0$, we can directly write such term for $\omega=-{\rm i}\kappa\neq-{\rm i}$ explicitly as a sum over residues, by deforming the contour around the positive real axis,
\begin{equation}
     \frac{1}{2\pi{\rm i}}\int_{\mathbb R} \frac{{\rm e}^{2{\rm i} B\omega'}G_-(\omega')}{G_+(\omega')d_+(\omega')d_-(\omega')}\frac{f(\omega)}{\omega'-{\rm i}\kappa} {\rm d} \omega' = -\sum_{n\ge 1} \frac{{\rm e}^{-2B\xi_n} {\rm i}\tilde\sigma_n f({\rm i}\xi_n)}{\xi_n-\kappa} .
\end{equation}
Plugging this result back into \eqref{epsilon_minus}, we find
\begin{multline}
    \frac{\epsilon_+({\rm i}\kappa)}{G_+({\rm i}\kappa)}=
    \frac{1}{2\pi{\rm i}}\int_{{\cal C}_\pm}\frac{{\rm e}^{-2B\xi}\delta\sigma({\rm i}\xi)v({\rm i}\xi)}{\xi-\kappa} {\rm d}\xi
    -\frac{m {\rm e}^ B}{2}\frac{G_+({\rm i})}{1+\kappa}\\
    +\frac{m{\rm e}^B}{2}\frac{1}{2\pi{\rm i}}\int_{{\cal C}_\pm} \frac{{\rm e}^{-2B\xi}G_+({\rm i}\xi)\delta\sigma({\rm i}\xi)}{(\xi-1)(\xi-\kappa)}{\rm d}\xi
    + \frac{h}{2\pi{\rm i}}\int_{{\cal C}_\pm}  \frac{G_+({\rm i}\xi)\delta\sigma({\rm i}\xi)(1-{\rm e}^{-2B\xi})}{\xi(\xi-\kappa)}{\rm d}\xi \\

   -\sum_{n\ge 1} \frac{{\rm e}^{-2B\xi_n} }{\xi_n-\kappa}({\rm i}\sigma_n^\pm-{\rm i}\tilde\sigma_n)\left(v({\rm i}\xi_n)-\frac{h G_+({\rm i}\xi_n)}{\xi_n}+\frac{m{\rm e}^BG_+({\rm i}\xi_n)}{2(\xi_n-1)}\right).
    \label{epsilon_plus_not1}
\end{multline}
We can see that in the boundary condition imposed by $\lim_{\kappa\rightarrow\infty}\kappa\epsilon_+({\rm i}\kappa)$ the only difference between the $\vartheta=0$ and $\vartheta=\pi$ cases are the residues, since ${\rm i}\tilde\sigma_n$ are absent in the former case. The remaining terms, which are identical in both, do not contain further non-perturbative corrections, other than those that come through $v$ from \eqref{vnotequal}.

As we saw in section \ref{sec_O3pi}, to calculate the free energy we need also to evaluate \eqref{epsilon_minus} at $\omega=-{\rm i}$. This  differs from \eqref{epsilon_plus_not1} because there is now an additional pole at $\omega={\rm i}$ when deforming some of the integrals. In the end we find
\begin{multline}
    \frac{\epsilon_+({\rm i})}{G_+({\rm i})}=

    \frac{1}{2\pi{\rm i}}\int_{{\cal C}_\pm} \frac{\delta\sigma({\rm i}\xi)}{\xi-1}\left({\rm e}^{-2B\xi}v({\rm i}\xi) +\frac{m{\rm e}^B{\rm e}^{-2B\xi}G_+({\rm i}\xi)}{2(\xi-1)}+\frac{h G_+({\rm i}\xi)(1-{\rm e}^{-2B\xi})}{\xi}\right) {\rm d}\xi
   \\
     -\frac{m {\rm e}^ B}{4}G_+({\rm i})
   -\sum_{n\ge 1} \frac{{\rm e}^{-2B\xi_n} ({\rm i}\sigma_n^\pm-{\rm i}\tilde\sigma_n)}{\xi_n-1}\left(v({\rm i}\xi_n)-\frac{h G_+({\rm i}\xi_n)}{\xi_n}+\frac{m{\rm e}^BG_+({\rm i}\xi_n)}{2(\xi_n-1)}\right)
   \\
   +\frac{m{\rm e}^{-B}}{2} \left({\rm i} G'_-({\rm i}\pm 0) -2 \Re\left\{{\rm i} G'_-({\rm i}) \right\}\right)
    \label{epsilon_1}.
\end{multline}
In the last term, we used \eqref{eq_d+d-_mod}. For $\vartheta=0$, as follows from the analysis of \cite{mmr-an}, we have
\begin{multline}
    \frac{\epsilon_+({\rm i})}{G_+({\rm i})}=

    \frac{1}{2\pi{\rm i}}\int_{{\cal C}_\pm} \frac{\delta\sigma({\rm i}\xi)}{\xi-1}
    \left({\rm e}^{-2B\xi}Q({\rm i}\xi) +\frac{m{\rm e}^B{\rm e}^{-2B\xi}G_+({\rm i}\xi)}{2(\xi-1)}+\frac{h G_+({\rm i}\xi)(1-{\rm e}^{-2B\xi})}{\xi}\right) {\rm d}\xi
   \\
   %
   -\frac{m {\rm e}^ B}{4}G_+({\rm i})
   -\! \sum_{n\ge 1} \frac{{\rm e}^{-2B\xi_n} {\rm i}\sigma_n^\pm}{\xi_n-1}
   \left( \! Q({\rm i}\xi_n)-\frac{h G_+({\rm i}\xi_n)}{\xi_n}+\frac{m{\rm e}^BG_+({\rm i}\xi_n)}{2(\xi_n-1)}\! \right)\!
   + \frac{m{\rm e}^{-B}}{2} {\rm i} G'_-({\rm i}\pm 0)
    \label{epsilon_0}.
\end{multline}
We can then see that the difference between these two equations is simply in the residues from the poles at ${\rm i}\xi_n$ as well as in the so-called IR pole. Since the first term in the second line is perturbative ($m{\rm e}^B \sim h$), the difference is then restricted to the non-perturbative terms, as expected.

In \eqref{epsilon_1}, it is already manifest that the change to the IR pole is that the real part changes sign while the ambiguous imaginary part remains the same. In addition, as a consequence of \eqref{eq_d+d-_mod}, we also have
\begin{equation}
  {\rm i}\sigma_n^\pm-{\rm i}\tilde\sigma_n = -\Re({\rm i}\sigma_n^\pm)+{\rm i}\Im({\rm i}\sigma_n^\pm),
  \label{res_ReIm}
\end{equation}
so the real part of the residues also changes sign.

Having found the map between the $\vartheta=0$ trans-series solution and that of $\vartheta=\pi$, we can apply it to our results in section \ref{sec_fen0}. The trans-series for the ground state energy is given by \eqref{egs_trans}, but with the constants
\begin{equation}
\begin{aligned}
{\cal C}_0^{\pm}(\pi) &=
\bigl[-\sin (\pi  \Delta )\mp{\rm i}  \cos (\pi  \Delta )\bigr]
\left(\frac{64}{{\rm e}^2}\right)^{\Delta }\frac{ \Gamma \left(\frac{1}{2}-\Delta \right) }{\left({\rm e} \Delta \right) \Gamma \left(\Delta +\frac{1}{2}\right)},\\
{\cal C}_1^{\pm}(\pi) &=
\left[\sin \biggl(\frac{\pi  \Delta }{1-2 \Delta }\biggr)\pm{\rm i} \cos \biggl(\frac{\pi  \Delta }{1-2 \Delta }\biggr)\right]
 \left(\frac{1024^{\Delta }}{4 {\rm e}^2 (1-2 \Delta )}\right)^{\frac{1}{1-2 \Delta }}
\frac{2 {\rm e} \left(1-2\Delta\right) \Gamma \bigl(\frac{1-4 \Delta }{2-4 \Delta }\bigr) }{ \Delta \Gamma \bigl(\frac{3-4 \Delta }{2-4 \Delta }\bigr)},
\end{aligned}
\label{fenCpm_pi}
\end{equation}
while the free energy can be written as
\begin{equation}
\label{fh-bos-pi}
\begin{aligned}
{\cal F} (h, \pi) &= - \frac{k^2 h^2}{4\tilde\alpha}\Biggl\{ 1
-
2\left(-\Re\left\{{\rm i} \sigma^+_1\right\}\pm{\rm i} \Im\left\{{\rm i} \sigma^+_1\right\}\right)\left({\rm e}^{-\frac{2}{\tilde{\alpha}}}\tilde{\alpha}^{1-2 \xi }\right)^{\xi _1}\!
\frac{\Bigl(\frac{G^2_+({\rm i})}{2 \pi  k^2}
\left(\frac{m}{\Lambda}\right)^{2}\Bigr)^{\xi_1}}{\left(\xi _1-1\right)^2 \xi _1}
+ \cdots \Biggr\}\\
&\qquad - \frac{m^2}{4\pi}  G_+({\rm i})\left(-\Re\left\{{\rm i} G'_-({\rm i})\right\}\pm{\rm i} \Im\left\{{\rm i} G'_-({\rm i}+ 0)\right\} \right).
\end{aligned}
\end{equation}
In either of these representations, one can see that the leading, real exponentially suppressed effects change sign, including $F(0, \pi)$, which now reads
\begin{equation}
F(0, \pi) = \frac{m^2}{16\sin(2\pi\Delta)}.
\end{equation}
Meanwhile, the formally imaginary terms, responsible for canceling the ambiguities induced by the large order behavior of the perturbative series,
remain unchanged.

In the $O(3)$ sigma model we identified the effects that change sign as due to instantons. In Fendley's coset sigma models
there are two types of corrections due to two types of poles: a sequence that survives the large $N$ limit $\Delta\rightarrow 0$, and a sequence that does not (this is similar to what happens in the PCF studied in \cite{mmr-an}).
Following the analysis of \cite{mmr-an}, we identify the first as ``renormalon-like'' and the second as ``instanton-like''. The preceding argument
about the change in the real part of the residues applies equally to both, i.e. the effect of the $\vartheta=\pi$ term on renormalon contributions is identical
to its effect on instanton contributions. Hence, our results provide evidence that,
in the presence of a non-trivial theta angle $\vartheta=\pi$, the real non-perturbative effects associated to renormalons change sign.

\subsection{Testing the analytic results}
\label{sec_test_results_Fendley}

As we did for the $O(3)$ sigma model, we now test the validity of the analytic expressions we obtained for the non-perturbative effects in
Fendley's integrable models. In particular, we check numerically the real exponential corrections to the normalized energy density $e/\rho^2$, as presented in \eqref{egs_trans}, with the coefficients of \eqref{fenCpm_0} (for $\vartheta=0$) and \eqref{fenCpm_pi} (for $\vartheta=\pi$). As in the case of the $O(3)$ sigma model, the numerical study of the integral equation \eqref{merTBA2} is subtle, and we used again the methods described in Appendix~\ref{app-numerical-integral}.

In \figref{fig_fullresultcompare_Fendley}, we plot the points $(\alpha,k^2 e/\rho^2)$ that we obtained numerically from the Bethe ansatz, by
using $\Delta=1/5$, at both $\vartheta=0$ and $\vartheta=\pi$. We compare the result to the Borel resummation of the perturbative series, shown as the dashed line. When $\vartheta=\pi$, we used the values $B=10/k$, with integers $1 \le k \le 10$. When $\vartheta=0$, we looked for the values of $B$ that
would result into the same numerical $\alpha$ that we found for $\vartheta=\pi$. To obtain the Borel resummation of the perturbative series, we used 81 coefficients, computed with the generalization of Volin's method.
\begin{figure}
\centering
\includegraphics[width=0.59\textwidth]{figures/fullresultcompare_Fendley_N5.pdf}
\caption{Numerical evaluation of $(\alpha,k^2 e/\rho^2\big)$ for Fendley's integrable model with $\Delta=1/5$, at $\vartheta=0$ (black dots) and $\vartheta=\pi$ (red dots). The dashed line corresponds to the Borel resummation of the perturbative part, $\alpha\Re\big[ s_{\pm}(\varphi)(\alpha) \big]$, which is shared by both values of the angle $\vartheta$.}
\label{fig_fullresultcompare_Fendley}
\end{figure}

In \figref{fig_e_rho_exp_Fendley} and Table \ref{tab_e_rho_exp_Fendley}, we consider the difference between the numerical result for
the Bethe ansatz, and the Borel resummation of the perturbative series.  We compare the result against the analytic prediction of \eqref{egs_trans},
and we find good agreement for $\Delta= 1/3$, $1/5$, $1/6$, at both $\vartheta=0$ and $\vartheta=\pi$. The leading exponential effect due to the IR pole,
represented by the dashed line, changes sign. This confirms the change of sign of the $F(0)$ term, i.e. the $h$-independent
term in the free energy. As for the subleading exponential effect, for $\Delta=1/5$ and $1/6$ it corresponds to a renormalon
of the type first identified in \cite{mmr-an}. The numerical test confirms that the real part of its contribution
does indeed change sign as we go from $\vartheta=0$ to $\vartheta=\pi$, as can be seen from the solid lines.
On the other hand, for $\Delta=1/3$, the singularity at $\xi_1$ corresponds to an instanton, and we confirm that
the associated contribution also changes sign, as expected. In this case, we mention that instead of the coefficient ${\cal C}_1^\pm$ in \eqref{fenCpm_0} and \eqref{fenCpm_pi}, which is only valid for the renormalon singularities, we had to use the the more general coefficient in \eqref{C1generic} with $\xi_1=3$ and its corresponding $\vartheta=\pi$ counterpart.%
\begin{figure}
\centering
\begin{subfigure}{\textwidth}
\includegraphics[trim= 0 0 1mm 0,width=0.43\textwidth]{figures/e_rho_exp_Fendley_N3.pdf}
\hspace{0.065\textwidth}
\includegraphics[width=0.43\textwidth]{figures/e_rho_exp_norm_Fendley_N3.pdf}
\caption{$\Delta=1/3$}
\vspace{0.75cm}
\end{subfigure}
\begin{subfigure}{\textwidth}
\includegraphics[trim= 0 0 1mm 0,width=0.43\textwidth]{figures/e_rho_exp_Fendley_N5.pdf}
\hspace{0.065\textwidth}
\includegraphics[width=0.43\textwidth]{figures/e_rho_exp_norm_Fendley_N5.pdf}
\caption{$\Delta=1/5$}
\vspace{0.75cm}
\end{subfigure}
\begin{subfigure}{\textwidth}
\includegraphics[trim= 0 0 1mm 0,width=0.43\textwidth]{figures/e_rho_exp_Fendley_N6.pdf}
\hspace{0.065\textwidth}
\includegraphics[width=0.43\textwidth]{figures/e_rho_exp_norm_Fendley_N6.pdf}
\caption{$\Delta=1/6$}
\end{subfigure}
\caption{In the figures on the left we plot the difference between the normalized energy density $k^2 e/\rho^2$ and the real part of the Borel resummation of the perturbative series, for the Fendley's integrable models at $\vartheta=0$ (black) and $\vartheta=\pi$ (red). The dots are the numerical
calculations of the difference. The line is the theoretical prediction, given by the exponential corrections of \eqref{egs_trans}, including the IR renormalon correction of order ${\rm e}^{-2/\alpha}$ (dashed line) and also the next subleading correction, of order ${\rm e}^{-2\xi_1/\alpha}$ (full line). In the figures on the right, the result has been normalized by dividing by the IR renormalon correction.}
\label{fig_e_rho_exp_Fendley}
\end{figure}

\begin{table}
\begin{equation*}
\arraycolsep=2.5mm
\begin{array}{llll}
\toprule
 \alpha & \text{Numerical } (\vartheta=0) & \text{Numerical } (\vartheta=\pi) &\text{Asymptotics}\\
\midrule
 0.090553 & (0\! \pm \! 6)\times 10^{-7} & (3\!\pm\! 18)\times 10^{-6} & \pm 2.32356 \times 10^{-10}\\
 0.169263 & 9.75(28\!\pm\! 32)\times 10^{-6} & -0.0000(05\!\pm\! 35) & \pm 9.75205 \times 10^{-6}\\
 0.240275 & 0.000392(06\! \pm \! 28 ) & -0.0003(8\!\pm\! 5) & \pm 0.000392\\
 0.305661 & 0.00263(31\! \pm \! 31) & -0.0026(3\!\pm\! 7) & \pm 0.002635\\
 0.366721 & 0.0084(80\! \pm \! 14) & -0.0084(8\!\pm\! 9) & \pm 0.008496\\
 0.424307 & 0.0187(7\!\pm\! 4) & -0.018(79\!\pm\! 11) & \pm 0.018811\\
 0.478970 & 0.0333(8\!\pm\! 9) & -0.033(41\!\pm\! 15) & \pm 0.033460\\
 0.531058 & 0.051(69\!\pm\! 15) & -0.051(65\!\pm\! 20) & \pm 0.051775\\
 0.580792 & 0.072(90\!\pm\! 24) & -0.072(64\!\pm\! 28) & \pm 0.072905\\
 0.628315 & 0.096(19\!\pm\! 34) & -0.095(5\!\pm\! 4) & \pm 0.096016\\
\bottomrule
\end{array}
\end{equation*}
\caption{Difference between the normalized energy density $k^2 e/\rho^2$ and the real
part of the Borel resummation of the perturbative series, for the Fendley's integrable model with $\Delta=1/5$. We compare the numerical value obtained from the Bethe ansatz (2nd and 3rd column), to the exponential corrections of \eqref{egs_trans} (4th column). The error estimates of $k^2 e/\rho^2$ are found by varying the quadrature parameters, see Appendix \ref{app-numerical-integral} for the details of the $\vartheta=\pi$ case.}
\label{tab_e_rho_exp_Fendley}
\end{table}

\sectiono{Discussion and prospects for future work}

\label{conclu}

In this paper we continued our analytic study of non-perturbative corrections in integrable field theories,
which we started in \cite{mmr-an}.
We considered integrable sigma models that admit a theta angle, and we studied the effect of this angle on
exponentially small corrections due to renormalons and instantons. As anticipated in the study of the deformed
sigma models in \cite{foz,fendley}, we found that would-be
instanton corrections change sign as we turn on $\vartheta= \pi$. We also found that, in some cases, the behavior of
renormalon corrections is more subtle: the imaginary parts of the residues associated to renormalon
singularities do not change sign, but their real parts do\footnote{In reaching this conclusion, we have
assumed a renormalon-like nature for these contributions, based for example on the fact that they survive the large $N$ limit.}.
This leads to renormalon contributions which are schematically of the form,
\begin{equation}
({\cal C}_1 (\vartheta)+ {\rm i} \, {\cal C}_2) g^b {\rm e}^{-A/g^2} \Phi(g),
\end{equation}
where $\Phi(g)$ is a formal power series and $b$ is an appropriate exponent. Most of the
ingredients of this trans-series, i.e. $A$, $b$ and $\Phi(g)$, can be in principle obtained from the large order
behavior of the perturbative series, through its imaginary part, proportional to ${\cal C}_2$. However, the trans-series parameter ${\cal C}_1 (\vartheta)$,
which is sensitive to the theta angle, cannot be predicted from perturbation theory alone. Such a scenario for renormalon
corrections was anticipated in e.g. \cite{grunberg}, and the models discussed in this paper
seem to provide a concrete realization thereof. The
possibility that the overall coefficients associated to renormalon corrections might depend on the theta angle, as we are discovering here,
was already pointed out in \cite{david-cpn}. More work is needed however to provide a microscopic description of the
non-perturbative corrections unveiled in these coset sigma models, and to verify the instanton/renormalon nature of the non-perturbative corrections.

Our analysis also establishes that purely real corrections which change sign as we change the theta angle can not
be detected in any way by a resurgent analysis of the perturbative series. These include notably the would-be instanton
correction of the $O(3)$ sigma model already discussed in \cite{mmr-an, bbh}.
We conclude that, in the terminology of \cite{dpmss, mmr-an}, these theories satisfy the weak version of the resurgence program
(physical observables are given by Borel-resummed trans-series), but not the strong version, since not all the trans-series corrections can be decoded from the perturbative series alone. This is however expected in theories which admit a theta angle (see \cite{dunne-unsal-cpn} for a similar discussion).

As we explained in the introduction, one important motivation for this work was to obtain exact results on non-perturbative effects
in semi-realistic asymptotically free theories, which can serve as signposts for non-perturbative methods in quantum field theory.
In the case of renormalon corrections, there is no known method to independently test the predictions we have obtained. However, the
analytic result (\ref{fh-ts}) provides a non-trivial prediction for an instanton correction in the $O(3)$ sigma model, which
could and should be testable by conventional methods. As we pointed out in (\ref{chith}), this correction leads to a non-vanishing topological
susceptibility in the presence of the $h$ field, and this prediction from integrability might be testable in a lattice calculation, along the lines of \cite{bnpw}.
Instanton calculus in non-supersymmetric theories
has been an ``arena for bloody controversies", in the famous description by Sidney Coleman \cite{coleman-as}, and we hope that this clean
prediction from integrability might help clarify some of the issues at stake. The presence of an external field $h$ should act as a
regulator for the various IR and UV singularities that afflict instanton calculations in the $O(3)$ sigma model, but understanding in a precise way
how this works is not obvious, at least for us.

\section*{Acknowledgements}
We would like to thank Ivan Graham and specially Martin Gander and Walter Gautschi for very useful
advice on the numerical study of the integral equation (\ref{merTBA2}). The numerical computations were performed at University of Geneva on ``Baobab" and ``Yggdrasil" HPC clusters. This work has been supported by the ERC-SyG project
``Recursive and Exact New Quantum Theory" (ReNewQuantum), which
received funding from the European Research Council (ERC) under the European
Union's Horizon 2020 research and innovation program,
grant agreement No. 810573.