\documentclass[
  aps,
  pre,
  preprint,
  longbibliography,
  floatfix
]{revtex4-2}

\usepackage[utf8]{inputenc}
\usepackage[T1]{fontenc}
\usepackage{newtxtext, newtxmath}
\usepackage[
  colorlinks=true,
  urlcolor=purple,
  citecolor=purple,
  filecolor=purple,
  linkcolor=purple
]{hyperref}
\usepackage{amsmath}
\usepackage{graphicx}
\usepackage{xcolor}
\usepackage{tikz-cd}
\usepackage[subfolder]{gnuplottex} % need to compile separately for APS
\usepackage{setspace}
\usepackage{tabularx}

\begin{document}

\title{Smooth and global Ising universal scaling functions}

\author{Jaron Kent-Dobias}
\affiliation{Laboratoire de Physique de l'Ecole Normale Supérieure, Paris, France}

\author{James P.~Sethna}
\affiliation{Laboratory of Atomic and Solid State Physics, Cornell University, Ithaca, NY, USA}

\date\today

\begin{abstract}
\end{abstract}

\maketitle

\section{Introduction}

At continuous phase transitions the thermodynamic properties of physical
systems have singularities. Celebrated renormalization group analyses imply
that not only the principal divergence but entire functions are
\emph{universal}, meaning that they will appear at any critical points that
connect phases of the same symmetries in the same spatial dimension. The study
of these universal functions is therefore doubly fruitful: it provides both a
description of the physical or model system at hand, and \emph{every other
system} whose symmetries, interaction range, and dimension puts it in the same
universality class.

The continuous phase transition in the two-dimensional Ising model is the most
well studied, and its universal thermodynamic functions have likewise received
the most attention. Precision numeric work both on lattice models and on the
``Ising'' conformal field theory (related by universality) have yielded
high-order polynomial expansions of those functions, along with a comprehensive
understanding of their analytic properties \cite{Fonseca_2003_Ising,
Mangazeev_2008_Variational, Mangazeev_2010_Scaling}.  In parallel, smooth
approximations of the Ising equation of state produce convenient, evaluable,
differentiable empirical functions \cite{Caselle_2001_The}. Despite being
differentiable, these approximations become increasingly poor when derivatives
are taken due to the absence of subtle singularities.

This paper attempts to find the best of both worlds: a smooth approximate
universal thermodynamic function that respects the global analyticity of the
Ising free energy. By constructing approximate functions with the correct
singularities, corrections converge exponentially to the true function. To make
the construction, we review the analytic properties of the Ising scaling
function. Parametric coordinates are introduced that remove unnecessary
singularities that are a remnant of the coordinates. Then, the arbitrary
analytic functions that compose those coordinates are approximated by truncated
polynomials whose coefficients are fixed by matching the series expansions of
the universal function.

\section{Universal scaling functions}

A renormalization group analysis predicts that certain thermodynamic functions
will be universal in the vicinity of any critical point in the Ising
universality class, from perturbed conformal fields to the end of the
liquid--gas coexistence line. Here we will review precisely what is meant by
universal.

Suppose one controls a temperature-like parameter $T$ and a magnetic field-like
parameter $H$, which in the proximity of a critical point at $T=T_c$ and $H=0$
have normalized reduced forms $t=(T-T_c)/T_c$ and $h=H/T$. Thermodynamic
functions are derived from the free energy per site $f=(F-F_c)/L^D$, which
depends on $t$, $h$, and a litany of irrelevant parameters we will henceforth
neglect.  Explicit renormalization with techniques like the
$\epsilon$-expansion or exact solutions like Onsager's can be used calculated
the flow of these parameters under continuous changes of scale $L=e^\ell$,
yielding equations of the form
\begin{align} \label{eq:raw.flow}
  \frac{dt}{d\ell}=\frac1\nu t+\cdots
  &&
  \frac{dh}{d\ell}=\frac{\beta\delta}\nu h+\cdots
  &&
  \frac{df}{d\ell}=Df+\cdots
\end{align}
where $D=2$ is the dimension of space and $\nu=1$, $\beta=\frac18$, and
$\delta=15$ are dimensionless constants. The combination
$\Delta=\beta\delta=\frac{15}8$ will appear often. The flow equations are
truncated here, but in general all terms allowed by the symmetries of the
parameters are present on their righthand side. By making a near-identity
transformation to the coordinates and the free energy of the form $u_t(t,
h)=t+\cdots$, $u_h(t, h)=h+\cdots$, and $u_f(f,t,h)=f+\cdots$, one can bring
the flow equations into the agreed upon simplest normal form
\begin{align} \label{eq:flow}
  \frac{du_t}{d\ell}=\frac1\nu u_t
  &&
  \frac{du_h}{d\ell}=\frac{\Delta}\nu u_h
  &&
  \frac{du_f}{d\ell}=Du_f-\frac1{4\pi}u_t^2
\end{align}
which are exact as written \cite{Raju_2019_Normal}. The flow of the
\emph{scaling fields} $u_t$ and $u_h$ is made exactly linear, while that of the
free energy is linearized as nearly as possible. The quadratic term in that
equation is unremovable due to a resonance between the value of $\nu$ and the
spatial dimension in two dimensions, while its coefficient is chosen as a
matter of convention, fixing the scale of $u_t$. Solving these equations for
$u_f$ yields
\begin{equation}
  \begin{aligned}
    u_f(u_t, u_h)
    &=|u_t|^{D\nu}\mathcal F_\pm(u_h|u_t|^{-\Delta})+\frac{|u_t|^{D\nu}}{8\pi}\log u_t^2 \\
    &=|u_h|^{D\nu/\Delta}\mathcal F_0(u_t|u_h|^{-1/\Delta})+\frac{|u_t|^{D\nu}}{8\pi}\log u_h^{2/\Delta} \\
  \end{aligned}
\end{equation}
where $\mathcal F_\pm$ and $\mathcal F_0$ are undetermined scaling functions.
The scaling functions are universal in the sense that if another system whose
critical point belongs to the same universality class has its parameters
brought to the form \eqref{eq:flow}, one will see the same functional form (up
to constant rescaling of $u_h$). The invariant scaling combinations that appear
as the arguments to the universal scaling functions will come up often, and we
will use $\xi=u_h|u_t|^{-\Delta}$ and $\eta=u_t|u_h|^{-1/\Delta}$.

The analyticity of the free energy at places away from the critical point implies that the functions
$\mathcal F_\pm$ and $\mathcal F_0$ have power-law expansions of their
arguments about zero. For instance, when $u_t$ goes to zero for nonzero $u_h$
there is no phase transition, and the free energy must be an analytic function
of its arguments. It follows that $\mathcal F_0$ is analytic about zero. This
is not the case at infinity: since $\mathcal F_0(\eta)=\eta^{D\nu}\mathcal
F_\pm(\eta^{-1/\Delta})$ has a power-law expansion about zero, $\mathcal
F_\pm(\xi)\sim \xi^{D\nu/\Delta}$ for large $\xi$. The nonanalyticity of
these functions at infinite argument can therefore be understood as an artifact
of the chosen coordinates.

For the scale of $u_t$ and $u_h$, we adopt the same convention as used by
\cite{Fonseca_2003_Ising}. The dependence of the nonlinear scaling variables on
the parameters $t$ and $h$ is system-dependent, and their form can be found for
common model systems (the square- and triangular-lattice Ising models) in the
literature \cite{Mangazeev_2010_Scaling, Clement_2019_Respect}.  To connect the
results of thes paper with Mangazeev and Fonseca, one can write $\mathcal
F_0(\eta)=\tilde\Phi(-\eta)=\Phi(-\eta)+(\eta^2/8\pi) \log \eta^2$ and
$\mathcal F_\pm(\xi)=G_{\mathrm{high}/\mathrm{low}}(\xi)$.

\section{Singularities}

\subsection{Essential singularity at the abrupt transition}

In the low temperature phase, the free energy has an essential singularity at
zero field, which becomes a branch cut along the negative-$h$ axis when
analytically continued to negative $h$ \cite{Langer_1967_Theory}. The origin
can be schematically understood to arise from a singularity that exists in the
imaginary free energy of the metastable phase of the model. When the
equilibrium Ising model with positive magnetization is subjected to a small
negative magnetic field, its equilibrium state instantly becomes one with a
negative magnetization. However, under physical dynamics it takes time to
arrive at this state, which happens after a fluctuation containing a
sufficiently large equilibrium `bubble' occurs.

The bulk of such a bubble of radius $R$ lowers the free energy by $2M|H|\pi
R^2$, where  $M$ is the magnetization, but its surface raises the free energy
by $2\pi R\sigma$, where $\sigma$ is the surface tension between the
stable--metastable interface. The bubble is sufficiently large to decay
metastable state when the differential bulk savings outweigh the surface costs.
This critical bubble occurs with free energy cost
\begin{equation}
  \Delta F_c
    \simeq\frac{\pi\sigma^2}{2M|H|}
    \simeq T\left(\frac{2M_0}{\pi\sigma_0^2}|\xi|\right)^{-1}
\end{equation}
where $\sigma_0$ and $M_0$ are the critical amplitudes for the surface tension
and magnetization, respectively \cite{Kent-Dobias_2020_Novel}.  In the context
of statistical mechanics, Langer demonstrated that the decay rate is
asymptotically proportional to the imaginary part of the free energy in the
metastable phase, with
\begin{equation}
  \operatorname{Im}F\propto\Gamma\sim e^{-\beta\Delta F_c}\simeq e^{-1/b|\xi|}
\end{equation}
which can be more rigorously related in the context of quantum field theory
[ref?]. The constant $b=2M_0/\pi\sigma_0^2$ is predicted by known properties,
and for our conventions for $u_t$ and $u_h$, $M_0/\sigma_0^2=\bar
s=2^{1/12}e^{-1/8}A^{3/2}$, where $A$ is Glaisher's constant.


\begin{figure}
  \includegraphics{figs/F_lower_singularities.pdf}
  \caption{
    Analytic structure of the low-temperature scaling function $\mathcal F_-$
    in the complex $\xi=u_h|u_t|^{-\Delta}\propto H$ plane. The circle
    depicts the essential singularity at the first order transition, while the
    solid line depicts Langer's branch cut.
  } \label{fig:lower.singularities}
\end{figure}

To lowest order, this singularity is a function of the scaling invariant $\xi$
alone. It is therefore suggestive that this should be considered a part of the
singular free energy and moreover part of the scaling function that composes
it. We will therefore make the ansatz that
\begin{equation} \label{eq:essential.singularity}
  \operatorname{Im}\mathcal F_-(\xi+i0)=A_0\Theta(-\xi)\xi e^{-1/b|\xi|}\left[1+O(\xi)\right]
\end{equation}
\cite{Houghton_1980_The}
The linear prefactor can be found through a more careful accounting of the
entropy of long-wavelength fluctuations in the droplet surface
\cite{Gunther_1980_Goldstone}. In conformal field theory, the prefactor is also known to be $A_0=\bar s/2\pi$.

\subsection{Yang--Lee edge singularity}

At finite size, the Ising model free energy is an analytic function of
temperature and field because it is the logarithm of a sum of positive analytic
functions. However, it can and does have singularities in the complex plane due
to zeros of the partition function at complex argument, and in particular at
imaginary values of field, $h$. Yang and Lee showed that in the thermodynamic
limit of the high temperature phase of the model, these zeros form a branch cut
along the imaginary $h$ axis that extends to $\pm i\infty$ starting at the
point $\pm ih_{\mathrm{YL}}$ \cite{Yang_1952_Statistical, Lee_1952_Statistical}.
The singularity of the phase transition occurs because these branch cuts
descend and touch the real axis as $T$ approaches $T_c$, with
$h_{\mathrm{YL}}\propto t^{\Delta}$. This implies that the
high-temperature scaling function for the Ising model should have complex
branch cuts beginning at $\pm i\xi_{\mathrm{YL}}$ for a universal constant
$\xi_{\mathrm{YL}}$.

\begin{figure}
  \includegraphics{figs/F_higher_singularities.pdf}
  \caption{
    Analytic structure of the high-temperature scaling function $\mathcal F_+$
    in the complex $\xi=u_h|u_t|^{-\Delta}\propto H$ plane. The squares
    depict the Yang--Lee edge singularities, while the solid lines depict
    branch cuts.
  } \label{fig:higher.singularities}
\end{figure}

The Yang--Lee singularities are critical points in their own right, with their
own universality class different from that of the Ising model
\cite{Fisher_1978_Yang-Lee}. Asymptotically close to this point, the scaling
function $\mathcal F_+$ takes the form
\begin{equation} \label{eq:yang.lee.sing}
  \mathcal F_+(\xi)
  =A(\xi) +B(\xi)[1+(\xi/\xi_{\mathrm{YL}})^2]^{1+\sigma}+\cdots
\end{equation}
with edge exponent $\sigma=-\frac16$ and $A$ and $B$ analytic functions at
$\xi_\mathrm{YL}$ \cite{Cardy_1985_Conformal, Fonseca_2003_Ising}. This creates
a branch cut stemming from the critical point along the imaginary-$\xi$ axis
with a growing imaginary part
\begin{equation}
  \operatorname{Im}\mathcal F_+(i\xi\pm0)=\pm A_\mathrm{YL}\frac12\Theta(\xi^2-\xi_\mathrm{YL}^2)[(\xi/\xi_\mathrm{YL})^2-1]^{1+\sigma}[1+O[(\xi-\xi_\mathrm{YL})^2]]
\end{equation}
This results in analytic structure for $\mathcal F_+$ shown in
Fig.~\ref{fig:higher.singularities}.

\section{Parametric coordinates}

The invariant combinations $\xi=u_h|u_t|^{-\Delta}$ or
$\eta=u_t|u_h|^{-1/\Delta}$ are natural variables to describe the scaling
functions, but prove unwieldy when attempting to make smooth approximations.
This is because, when defined in terms of these variables, scaling functions
that have polynomial expansions at small argument have nonpolynomial expansions
at large argument. Rather than deal with the creative challenge of dreaming up
functions with different asymptotic expansions in different limits, we adopt
another coordinate system, in terms of which a scaling function can be defined
that has polynomial expansions in \emph{all} limits.

In all dimensions, the Schofield coordinates $R$ and $\theta$ will be implicitly defined by
\begin{align} \label{eq:schofield}
  u_t(R, \theta) = R(1-\theta^2)
  &&
  u_h(R, \theta) = R^{\Delta}g(\theta)
\end{align}
where $g$ is an odd function whose first zero lies at $\theta_0>1$. We take
\begin{align} \label{eq:schofield.funcs}
  g(\theta)=\left(1-\frac{\theta^2}{\theta_0^2}\right)\sum_{i=0}^\infty g_i\theta^{2i+1}.
\end{align}
This means that $\theta=0$ corresponds to the high-temperature zero-field line,
$\theta=1$ to the critical isotherm at nonzero field, and $\theta=\theta_0$ to
the low-temperature zero-field (phase coexistence) line.
In practice the infinite series in \eqref{eq:schofield.funcs} cannot be
entirely fixed, and it will be truncated at finite order.

One can now see the convenience of these coordinates. Both invariant scaling
combinations depend only on $\theta$, as
\begin{align}
  \xi=u_h|u_t|^{-\Delta}=\frac{g(\theta)}{|1-\theta^2|^{\Delta}} &&
  \eta=u_t|u_h|^{-1/\Delta}=\frac{1-\theta^2}{|g(\theta)|^{1/\Delta}}
\end{align}
Moreover, both scaling variables have polynomial expansions in $\theta$ near zero, with
\begin{align}
  &\xi= g'(0)\theta+\cdots  && \text{for $\theta\simeq0$}\\
  &\xi=g'(\theta_0)(\theta_0^2-1)^{-\Delta}(\theta-\theta_0)+\cdots && \text{for $\theta\simeq\theta_0$}
  \\
  &\eta=-2(\theta-1)g(1)^{-1/\Delta}+\cdots && \text{for $\theta\simeq1$}.
\end{align}
Since the scaling functions $\mathcal F_\pm(\xi)$ and $\mathcal F_0(\eta)$ have
polynomial expansions about small $\xi$ and $\eta$, respectively, this implies
both will have polynomial expansions in $\theta$ everywhere.

Therefore, in Schofield coordinates one expects to be able to define a global
scaling function $\mathcal F(\theta)$ which has a polynomial expansion in its
argument for all real $\theta$ by
\begin{equation}
  u_f(R,\theta)=R^{D\nu}\mathcal F(\theta)+(1-\theta^2)^2\frac{R^2}{8\pi}\log R^2
\end{equation}
For small $\theta$, $\mathcal F(\theta)$ will
resemble $\mathcal F_+$, for $\theta$ near one it will resemble $\mathcal F_0$,
and for $\theta$ near $\theta_0$ it will resemble $\mathcal F_-$. This can be seen explicitly using the definitions \eqref{eq:schofield} to relate the above form to the original scaling functions, giving
\begin{equation} \label{eq:scaling.function.equivalences.2d}
  \begin{aligned}
    \mathcal F(\theta)
    &=|t(\theta)|^{D\nu}\mathcal F_\pm\left[g(\theta)|1-\theta^2|^{-\Delta}\right]
    +\frac{(1-\theta^2)^2}{8\pi}\log t(\theta)^2\\
    &=|h(\theta)|^{D\nu/\Delta}\mathcal F_0\left[(1-\theta^2)|g(\theta)|^{-1/\Delta}\right]
    +\frac{(1-\theta^2)^2}{8\pi}\log g(\theta)^{2/\Delta}
  \end{aligned}
\end{equation}
This leads us
to expect that the singularities present in these functions will likewise be
present in $\mathcal F(\theta)$. The analytic structure of this function is
shown in Fig.~\ref{fig:schofield.singularities}. Two copies of the Langer
branch cut stretch out from $\pm\theta_0$, where the equilibrium phase ends,
and the Yang--Lee edge singularities are present on the imaginary-$\theta$
line, where they must be since $\mathcal F$ has the same symmetry in $\theta$
as $\mathcal F_+$ has in $\xi$.

The location of the Yang--Lee edge singularities can be calculated directly
from the coordinate transformation \eqref{eq:schofield}. Since $g(\theta)$ is
an odd real polynomial for real $\theta$, it is imaginary for imaginary
$\theta$. Therefore,
\begin{equation}
  i\xi_{\mathrm{YL}}=\frac{g(i\theta_{\mathrm{YL}})}{(1+\theta_{\mathrm{YL}}^2)^{-\Delta}}
\end{equation}
The location $\theta_0$ is not fixed by any principle.

\begin{figure}
  \includegraphics{figs/F_theta_singularities.pdf}
  \caption{
    Analytic structure of the global scaling function $\mathcal F$ in the
    complex $\theta$ plane. The circles depict essential singularities of the
    first order transitions, the squares the Yang--Lee singularities, and the
    solid lines depict branch cuts.
  } \label{fig:schofield.singularities}
\end{figure}

\section{Functional form for the parametric free energy}

As we have seen in the previous sections, the unavoidable singularities in the
scaling functions are readily expressed as singular functions in the imaginary
part of the free energy.

Our strategy follows. First, we take the singular imaginary parts of the
scaling functions $\mathcal F_{\pm}(\xi)$ and truncate them to the lowest order
accessible under polynomial coordinate changes of $\xi$. Then, we assert that
the imaginary part of $\mathcal F(\theta)$ must have this simplest form,
implicitly defining the parametric coordinate change. Third, we perform a
Kramers--Kronig type transformation to establish an explicit form for the real
part of $\mathcal F(\theta)$. Finally, we make good on the assertion made in
the second step by finding the coordinate transformation that produces the
correct series coefficients of $\mathcal F_{\pm}$.

This success of this stems from the commutative diagram below. So long as the
application of Schofield coordinates and the Kramers--Kronig relation can be
said to commute, we may assume we have found correct coordinates for the
simplest form of the imaginary part to be fixed later by the real part.
\[
  \begin{tikzcd}[row sep=large, column sep = 9em]
  \operatorname{Im}\mathcal F_\pm(\xi) \arrow{r}{\text{Kramers--Kronig in $\xi$}} \arrow[]{d}{\text{Schofield}} & \operatorname{Re}\mathcal F_{\pm}(\xi) \arrow{d}{\text{Schofield}} \\%
  \operatorname{Im}\mathcal F(\theta) \arrow{r}{\text{Kramers--Kronig in $\theta$}}& \operatorname{Re}\mathcal F(\theta)
\end{tikzcd}
\]
We require that, for $\theta\in\mathbb R$
\begin{equation} \label{eq:imaginary.abrupt}
  \operatorname{Im}\mathcal F(\theta+0i)=\operatorname{Im}\mathcal F_0(\theta+0i)=C_0[\Theta(\theta-\theta_0)\mathcal I(\theta)-\Theta(-\theta-\theta_0)\mathcal I(-\theta)]
\end{equation}
where
\begin{equation}
  \mathcal I(\theta)=(\theta-\theta_0)e^{-1/B(\theta-\theta_0)}
\end{equation}
reproduces the essential singularity in \eqref{eq:essential.singularity}.
Independently, we require for $\theta\in\mathbb R$
\begin{equation}
  \operatorname{Im}\mathcal F(i\theta+0)=\operatorname{Im}\mathcal F_\mathrm{YL}(i\theta+0)=\frac{C_\mathrm{YL}}2\Theta(\theta^2-\theta_\mathrm{YL}^2)[(\theta/\theta_\mathrm{YL})^2-1]^{1+\sigma}
\end{equation}
Fixing these requirements for the imaginary part of $\mathcal F(\theta)$ fixes
its real part up to an analytic even function $G(\theta)$, real for real $\theta$.

\begin{figure}
  \includegraphics{figs/contour_path.pdf}
  \caption{
    Integration contour over the global scaling function $\mathcal F$ in the
    complex $\theta$ plane used to produce the dispersion relation. The
    circular arc is taken to infinity, while the circles around the
    singularities are taken to zero.
  } \label{fig:contour}
\end{figure}

To find the real part of the nonanalytic part of the scaling function, we make
use of the identity
\begin{equation}
  0=\oint_{\mathcal C}d\vartheta\,\frac{\mathcal F(\vartheta)}{\vartheta^2(\vartheta-\theta)}
\end{equation}
where $\mathcal C$ is the contour in Figure \ref{fig:contour}. The integral is
zero because there are no singularites enclosed by the contour. The only
nonvanishing contributions from this contour as the radius of the semicircle is
taken to infinity are along the real line and along the branch cut in the upper
half plane. For the latter contributions, the real parts of the integration up
and down cancel out, while the imaginary part doubles. This gives
\begin{equation}
  \begin{aligned}
    0&=\left[\int_{-\infty}^\infty+\lim_{\epsilon\to0}\left(\int_{i\infty-\epsilon}^{i\theta_{\mathrm{YL}}-\epsilon}+\int^{i\infty+\epsilon}_{i\theta_{\mathrm{YL}}+\epsilon}\right)\right]
      d\vartheta\,\frac{\mathcal F(\vartheta)}{\vartheta^2(\vartheta-\theta)} \\
     &=\int_{-\infty}^\infty d\vartheta\,\frac{\mathcal F(\vartheta)}{\vartheta^2(\vartheta-\theta)}
     +2i\int_{i\theta_{\mathrm{YL}}}^{i\infty}d\theta'\,\frac{\operatorname{Im}\mathcal F(\vartheta)}{\vartheta^2(\vartheta-\theta)} \\
     &=-i\pi\frac{\mathcal F(\theta)}{\theta^2}+\mathcal P\int_{-\infty}^\infty d\vartheta\,\frac{\mathcal F(\vartheta)}{\vartheta^2(\vartheta-\theta)}
     +2i\int_{i\theta_{\mathrm{YL}}}^{i\infty}d\vartheta\,\frac{\operatorname{Im}\mathcal F(\vartheta)}{\vartheta^2(\vartheta-\theta)}
  \end{aligned}
\end{equation}
where $\mathcal P$ is the principle value.  In principle one would need to
account for the residue of the pole at zero, but since its order is less than
two and $\mathcal F(0)=\mathcal F'(0)=0$, this evaluates to zero. Rearranging, this gives
\begin{equation}
  \mathcal F(\theta)
  =\frac{\theta^2}{i\pi}\mathcal P\int_{-\infty}^\infty d\vartheta\,\frac{\mathcal F(\vartheta)}{\vartheta^2(\vartheta-\theta)}
  +\frac{2\theta^2}\pi\int_{i\theta_{\mathrm{YL}}}^{i\infty}d\vartheta\,\frac{\operatorname{Im}\mathcal F(\theta')}{\vartheta^2(\vartheta-\theta)}
\end{equation}
Taking the real part of both sides, we find
\begin{equation}
  \operatorname{Re}\mathcal F(\theta)
  =\frac{\theta^2}{\pi}\mathcal P\int_{-\infty}^\infty d\vartheta\,\frac{\operatorname{Im}\mathcal F(\vartheta)}{\vartheta^2(\vartheta-\theta)}
  -\frac{2\theta^2}\pi\int_{\theta_{\mathrm{YL}}}^{\infty}d\vartheta\,\frac{\operatorname{Im}\mathcal F(i\vartheta)}{\vartheta(\vartheta^2+\theta^2)}
\end{equation}
Because the real part of $\mathcal F$ is even, the imaginary part must be odd. Therefore
\begin{equation} \label{eq:dispersion}
  \operatorname{Re}\mathcal F(\theta)
  =\frac{\theta^2}{\pi}
  \int_{\theta_0}^\infty d\vartheta\,\frac{\operatorname{Im}\mathcal F(\vartheta)}{\vartheta^2}\left(\frac1{\vartheta-\theta}+\frac1{\vartheta+\theta}\right)
  -\frac{2\theta^2}\pi\int_{\theta_{\mathrm{YL}}}^{\infty}d\vartheta\,\frac{\operatorname{Im}\mathcal F(i\vartheta)}{\vartheta(\vartheta^2+\theta^2)}
\end{equation}
Evaluating these ordinary integrals, we find for $\theta\in\mathbb R$
\begin{equation}
  \operatorname{Re}\mathcal F(\theta)=\operatorname{Re}\mathcal F_0(\theta)+\mathcal F_\mathrm{YL}(\theta)+G(\theta)
\end{equation}
where
\begin{equation} \label{eq:2d.real.Fc}
  \operatorname{Re}\mathcal F_0(\theta)
  =C_0[\mathcal R(\theta)+\mathcal R(-\theta)]
\end{equation}
where $\mathcal R$ is given by the function
\begin{equation}
  \mathcal R(\theta)
  =\frac1\pi\left[
    \theta_0e^{1/B\theta_0}\operatorname{Ei}(-1/B\theta_0)
    +(\theta-\theta_0)e^{-1/B(\theta-\theta_0)}\operatorname{Ei}(1/B(\theta-\theta_0))
  \right]
\end{equation}
and
\begin{equation}
  \mathcal F_{\mathrm{YL}}(\theta)=C_{\mathrm{YL}}\left[(\theta^2+\theta_{\mathrm{YL}}^2)^{1+\sigma}-\theta_{\mathrm{YL}}^{2(1+\sigma)}\right]
\end{equation}
We have also included the analytic part $G$, which we assume has a simple series expansion
\begin{equation}
  G(\theta)=\sum_{i=1}^\infty G_i\theta^{2i}
\end{equation}
From the form of the real part, we can infer the form of $\mathcal F$ that is analytic for the whole complex plane except at the singularities and branch cuts previously discussed.
For $\theta\in\mathbb C$, we take
\begin{equation}
  \mathcal F(\theta)=\mathcal F_0(\theta)+\mathcal F_{\mathrm{YL}}(\theta)+G(\theta),
\end{equation}
where
\begin{equation}
  \mathcal F_0(\theta)=C_0\left\{
    \mathcal R(\theta)+\mathcal R(-\theta)+i\operatorname{sgn}(\operatorname{Im}\theta)[\mathcal I(\theta)-\mathcal I(-\theta)]
  \right\}
\end{equation}

\section{Fitting}

The scaling function has a number of free parameters: the position $\theta_0$
of the abrupt transition, prefactors in front of singular functions from the
abrupt transition and the Yang--Lee point, the coefficients in the analytic
part $G$ of $\mathcal F$, and the coefficients in the undetermined function
$g$.  Other parameters are determined by known properties.

For $\theta>\theta_0$, the form \eqref{eq:essential.singularity} can be
expanded around $\theta=\theta_0$ to yield
\begin{equation}
  \begin{aligned}
    \operatorname{Im}u_f
    &\simeq A_0 u_t(\theta)^{D\nu}\xi(\theta)\exp\left\{\frac1{b\xi(\theta)}\right\} \\
    &=A_0R^{D\nu}(\theta_0^2-1)^{D\nu}\xi'(\theta_0)(\theta-\theta_0)
    \exp\left\{\frac1{b\xi'(\theta_0)}\left(\frac1{\theta-\theta_0}
      -\frac{\xi''(\theta_0)}{2\xi'(\theta_0)}\right)
      \right\}\left(1+O[(\theta-\theta_0)^2]\right)
  \end{aligned}
\end{equation}
Comparing this with the requirement \eqref{eq:imaginary.abrupt}, we find that
\begin{equation}
  B=-b\xi'(\theta_0)=-b\frac{g'(\theta_0)}{(\theta_0^2-1)^{1/\Delta}}
\end{equation}
and
\begin{equation}
  \begin{aligned}
    C_0&=A_0t(\theta_0^2-1)^{D\nu}\xi'(\theta_0)\exp\left\{
    -\frac{\xi''(\theta_0)}{2b\xi'(\theta_0)^2}
  \right\} \\
       &=
       A_0(\theta_0^2-1)^{D\nu-\Delta}g'(\theta_0)
       \exp\left\{-\frac1b\left(\frac{(\theta_0^2-1)^\Delta g''(\theta_0)}{2g'(\theta_0)^2}-\frac{2\Delta(\theta_0^2-1)^{\Delta - 1}\theta_0}{g'(\theta_0)}
       \right)\right\}
  \end{aligned}
\end{equation}
fixing $B$ and $C_0$.

This leaves as unknown variables the positions $\theta_0$ and
$\theta_{\mathrm{YL}}$ of the abrupt transition and Yang--Lee edge singularity,
the amplitude $C_\mathrm{YL}$ of the latter, and the unknown functions $G$ and
$g$. We determine these approximately by iteration in the polynomial order at
which the free energy and its derivative matches known results. We write as a
cost function the difference between the known series coefficients of the
scaling functions $\mathcal F_\pm$ and the series coefficients of our
parametric form evaluated at the same points, $\theta=0$ and $\theta=\theta_0$,
weighted by the uncertainty in the value of the known coefficients or by a
machine-precision cutoff, whichever is larger. A Levenburg--Marquardt algorithm
is performed on the cost function to find a parameter combination which
minimizes it. As larger polynomial order in the series are fit, the truncations
of $F$ and $h$ are extended to higher order so that the codimension of the fit
is constant. A term is added to $F$ whenever a new coefficient of the high
temperature series is added, and one is added to $h$ whenever a new coefficient
of the low temperature series is added.

We performed this procedure starting with $n=2$, or matching the scaling
function at the low and high temperature zero field points to quadratic order,
through $n=9$. The resulting fit coefficients can be found in Table
\ref{tab:fits} without any sort of uncertainty, which is difficult to quantify
directly due to the truncation of series. However, precise results exist for
the value of the scaling function at the critical isotherm, or equivalently for
the series coefficients of the scaling function $\mathcal F_0$, and the
accuracy of the fit results can be checked against the known values here.

\begin{table}\label{tab:fits}
  \begin{tabular}{c|ccc}
    $n$ & $\mathcal F_-^{(n)}$ & $\mathcal F_0^{(n)}$ & $\mathcal F_+^{(n)}$ \\\hline
    0   & 0                    & $-1.197733383797993$ & 0                    \\
    1   & $-1.35783834$        & $-0.318810124891$    & 0                    \\
    2   & $-0.048953289720$    & $0.110886196683$     & $-1.84522807823$     \\
    3   & 0.0388639290         & $0.01642689465$      & 0                    \\
    4   & $-0.068362121$       & $-2.639978\times10^{-4}$ & 8.3337117508     \\
    5   & 0.18388371           & $-5.140526\times10^{-4}$ & 0                \\
    6   & $-0.659170$          & $2.08856\times 10^{-4}$ & $-95.16897$       \\
    7   & 2.937665             & $-4.4819\times10^{-5}$  & 0                 \\
    8   & $-15.61$             & $3.16\times10^{-7}$  & 1457.62              \\
    9   & 96.76                & $4.31\times10^{-6}$  & 0                    \\
    10  & $-679$               & $-1.99\times10^{-6}$ & -25891               \\
    11  & $5.34\times10^3$     &                      & 0                    \\
    12  & $-4.66\times10^4$    &                      & $5.02\times10^5$     \\
    13  & $4.46\times10^5$     &                      & 0                    \\
    14  & $-4.66\times10^6$    &                      & $-1.04\times10^7$
  \end{tabular}
\end{table}

\begin{table}
  \singlespacing
  \begin{tabular}{c|llllllll}
    \multicolumn{1}{c|}{$n$} &
     \multicolumn{1}{c}{$\theta_\mathrm{YL}$} &
     \multicolumn{1}{c}{$A_\mathrm{YL}$} &
     \multicolumn{1}{c}{$F_1$} &
     \multicolumn{1}{c}{$F_2$} &
     \multicolumn{1}{c}{$F_3$} &
     \multicolumn{1}{c}{$F_4$} &
     \multicolumn{1}{c}{$F_5$} &
     \multicolumn{1}{c}{$F_6$} \\
    \hline
      2 &
      0.18041 &
      2.1295 &
      1.2447 &
      0.49975 \\
      3 &
      0.19613 &
      2.1999 &
      1.2402 &
      0.39028 \\
      4 &
      0.19563 &
      2.2433 &
      1.2964 &
      0.37080 &
      $-0.028926$ \\
      5 &
      0.19553 &
      2.2321 &
      1.2804 &
      0.36318 &
      $-0.028924$ & & \\
      6 &
      0.19737 &
      2.3981 &
      1.5091 &
      0.39200 &
      $-0.090023$ &
      0.017233 & \\
      7 &
      0.19730 &
      2.4229 &
      1.5454 &
      0.41320 &
      $-0.12161$ &
      0.026346 & \\
      8 &
      0.19655 &
      2.5513 &
      1.7323 &
      0.59677 &
      $-0.29521$ &
      0.078509 &
      $-0.0072514$ \\
      9 &
      1.3754 &
      0.19652 &
      2.5482 &
      1.7278 &
      0.60278 &
      $-0.28737$ &
      0.072411 &
      $-0.0072455$ \\
      \hline
    \end{tabular}
    \begin{tabular}{c|llllllll}
        \hline
      $n$ &
        \multicolumn{1}{c}{$\theta_0$} &
        \multicolumn{1}{c}{$h_1$} &
        \multicolumn{1}{c}{$h_2$} &
        \multicolumn{1}{c}{$h_3$} &
        \multicolumn{1}{c}{$h_4$} &
        \multicolumn{1}{c}{$h_5$} &
        \multicolumn{1}{c}{$h_6$} &
        \multicolumn{1}{c}{$h_7$} \\
      \hline
        2 &
        1.2114 \\
        3 &
      1.3498 &
      $-0.014909$ \\
      4 &
      1.4490 &
      $-0.10871$ & $-0.0031747$ \\
      5 &
      1.4719 &
      $-0.11399$ & $-0.0031669$ & $8.8574\times10^{-7}$ \\
      6 &
      1.4358 &
      $-0.19533$ & 0.029301 & 0.0039906 & $-0.00011913$ \\
      7 &
      1.4324 &
      $-0.22077$ & 0.036245 & 0.010120 & $-0.0011434$ & 0.00010095 \\
      8 &
      1.3710 &
      $-0.35150$ & 0.0050232 & 0.053659 & $-0.019806$ & 0.0033531 & $-0.00026034$ \\
  \end{tabular}
\end{table}

\begin{figure}
  \begin{gnuplot}[terminal=epslatex]
    dat = 'data/phi_comparison.dat'

    set xlabel '$n$'
    set ylabel '$|\Delta\mathcal F_0^{(m)}(0)|$'

    set format y '$10^{%T}$'
    set style data linespoints
    set logscale y
    set key title '\raisebox{0.5em}{$m$}' bottom left
    set xrange [1.5:9.5]

    plot \
      dat using 1:2 title '0', \
      dat using 1:3 title '1', \
      dat using 1:4 title '2', \
      dat using 1:5 title '3', \
      dat using 1:6 title '4'
  \end{gnuplot}
  \caption{
    The error in the $m$th derivative of the scaling function $\mathcal F_0$
    with respect to $\eta$ evaluated at $\eta=0$, as a function of the
    polynomial order $n$ at which the scaling function was fit.
  }
\end{figure}

\begin{figure}
  \begin{gnuplot}[terminal=epslatex]
    dat = 'data/yl_comparison.dat'

    set xlabel '$n$'
    set ylabel '$|\Delta\xi_\mathrm{YL}|$'
    set xrange [1.5:9.5]
    set yrange [0.000005:0.05]

    set format y '$10^{%T}$'
    set style data yerrorlines
    set logscale y
    unset key

    plot dat using 1:2:3
  \end{gnuplot}
  \caption{
    The error in the location of the Yang--Lee edge singularity as a function
    of the polynomial order $n$ at which the scaling function was fit. Error
    bars denote the uncertainty in the known location of the singularity.
  }
\end{figure}

\begin{figure}
  \begin{gnuplot}[terminal=epslatex, terminaloptions={size 8.65cm,5.35cm}]
    dat9 = 'data/h_series_ours_9.dat'
    dat11 = 'data/h_series_ours_11.dat'
    dat13 = 'data/h_series_ours_13.dat'
    dat15 = 'data/h_series_ours_15.dat'

    ratLast(x) = (back2 = back1, back1 = x, back1 / back2)

    back1 = 0
    back2 = 0

    set xrange [0:1.05]
    set yrange [0:0.55]

    set xlabel '$1/n$'
    set ylabel '$h_n/h_{n-1}$'

    set style data linespoints

    plot \
      dat9 using (1/($0)):(abs(ratLast($1))) title '9', \
      dat11 using (1/($0)):(abs(ratLast($1))) title '11', \
      dat13 using (1/($0)):(abs(ratLast($1))) title '13', \
      dat15 using (1/($0)):(abs(ratLast($1))) title '15', \
      0.5 - 0.675 * x lc black
  \end{gnuplot}
  \caption{
  }
\end{figure}

\subsection{Comparison}

\begin{figure}
  \begin{gnuplot}[terminal=epslatex]
    dat1 = 'data/glow_series_numeric.dat'
    dat2 = 'data/glow_series_ours_0.dat'
    dat3 = 'data/glow_series_caselle.dat'
    set key top left Left reverse
    set logscale y
    set xlabel '$n$'
    set ylabel '$\mathcal F_-^{(n)}$'
    set format y '$10^{%T}$'

    plot \
      dat1 using 1:(abs($2)) title 'Numeric' with yerrorbars, \
      dat2 using 1:(abs($2)) title 'Ours ($n=2$)', \
      dat3 using 1:(abs($2)) title 'Caselle'
  \end{gnuplot}
  \caption{
  }
\end{figure}

\begin{figure}
  \begin{gnuplot}[terminal=epslatex]
    dat1 = 'data/glow_series_numeric.dat'
    dat2 = 'data/glow_series_ours_0.dat'
    dat3 = 'data/glow_series_caselle.dat'
    ratLast(x) = (back2 = back1, back1 = x, back1 / back2)
    back1 = 0
    back2 = 0
    set xlabel '$1/n$'
    set xrange [0:0.55]
    set ylabel '$\mathcal F_n/\mathcal F_{n-1}$'
    set yrange [0:15]

    plot \
      dat1 using (1/$1):(abs(ratLast($2))) title 'Numeric', \
      dat2 using (1/$1):(abs(ratLast($2))) title 'Ours ($n=2$)', \
      dat3 using (1/$1):(abs(ratLast($2))) title 'Caselle'
  \end{gnuplot}
  \caption{
  }
\end{figure}

\begin{figure}
  \begin{gnuplot}[terminal=epslatex]
    dat1 = 'data/ghigh_numeric.dat'
    set key top left Left reverse
    set logscale y
    set xlabel '$n$'
    set ylabel '$\mathcal F_+^{(n)}$'
    set format y '$10^{%T}$'

    plot \
      dat1 using 1:(abs($2)) title 'Numeric' with yerrorbars
  \end{gnuplot}
  \caption{
  }
\end{figure}

\begin{figure}
  \begin{gnuplot}[terminal=epslatex]
    dat1 = 'data/phi_numeric.dat'
    set key top left Left reverse
    set logscale y
    set xlabel '$n$'
    set ylabel '$|\mathcal F_0^{(n)}|$'
    set format y '$10^{%T}$'
    set xrange [-0.5:10.5]

    plot \
      dat1 using 1:(abs($2)) title 'Numeric' with yerrorbars
  \end{gnuplot}
  \caption{
  }
\end{figure}

\section{Outlook}

The successful smooth description of the Ising free energy produced in part by analytically continuing the singular imaginary part of the metastable free energy inspires an extension of this work: a smooth function that captures the universal scaling \emph{through the coexistence line and into the metastable phase}. Indeed, the tools exist to produce this: by writing $t(\theta)=(1-\theta^2)(1-(\theta/\theta_m)^2)$ for some $\theta_m>\theta_0$, the invariant scaling combination

\begin{acknowledgments}
  The authors would like to thank Tom Lubensky, Andrea Liu, and Randy Kamien
  for helpful conversations. The authors would also like to think Jacques Perk
  for pointing us to several insightful studies. JPS thanks Jim Langer for past
  inspiration, guidance, and encouragement. This work was supported by NSF
  grants DMR-1312160 and DMR-1719490.
\end{acknowledgments}

\bibliography{ising_scaling}

\end{document}