%
%  Created by Jaron Kent-Dobias on Thu Apr 20 12:50:56 EDT 2017.
%  Copyright (c) 2017 Jaron Kent-Dobias. All rights reserved.
%
\documentclass[aps,prl,reprint]{revtex4-1}

\usepackage[utf8]{inputenc}
\usepackage{amsmath,amssymb,latexsym,mathtools,xifthen}

% uncomment to label only equations that are referenced in the text
%\mathtoolsset{showonlyrefs=true}

% I want labels but don't want to type out ``equation''
\def\[{\begin{equation}}
\def\]{\end{equation}}

% math not built-in
\def\arcsinh{\mathop{\mathrm{arcsinh}}\nolimits}
\def\arccosh{\mathop{\mathrm{arccosh}}\nolimits}
\def\ei{\mathop{\mathrm{Ei}}\nolimits} % exponential integral Ei
\def\re{\mathop{\mathrm{Re}}\nolimits}
\def\im{\mathop{\mathrm{Im}}\nolimits}
\def\sgn{\mathop{\mathrm{sgn}}\nolimits}
\def\dd{d} % differential
\def\O{O}          % big O
\def\o{o}          % little O

% subscript for ``critical'' values, e.g., T_\c
\def\c{\mathrm c}

% scaling functions
\def\fM{\mathcal M}  % magnetization
\def\fX{\mathcal Y}  % susceptibility
\def\fXt{\check{\mathcal Y}}  % susceptibility
\def\fF{\mathcal F}  % free energy
\def\fiF{\mathcal H} % imaginary free energy
\def\fS{\mathcal S}  % surface tension
\def\fG{\mathcal G}  % exponential factor

% lattice types
\def\sq{\mathrm{sq}}
\def\tri{\mathrm{tri}}
\def\hex{\mathrm{hex}}

% dimensions
\def\dim{d}
\def\twodee{\textsc{2d}}
\def\threedee{\textsc{3d}}
\def\fourdee{\textsc{4d}}

% fancy partial derivative
\newcommand\pd[3][]{
  \ifthenelse{\isempty{#1}}
    {\def\tmp{}}
    {\def\tmp{^#1}}
  \frac{\partial\tmp#2}{\partial#3\tmp}
}

% used to reformat display math to fit in two-column ``reprint'' mode
\makeatletter
\newif\ifreprint
\@ifclasswith{revtex4-1}{reprint}{\reprinttrue}{\reprintfalse}
\makeatother

\begin{document}

\title{Essential Singularities in Universal Scaling Functions at the Ising Coexistence Line}
\author{Jaron Kent-Dobias}
\author{James P.~Sethna}
\affiliation{Laboratory of Atomic and Solid State Physics, Cornell University, Ithaca, NY, USA}

\date\today

\begin{abstract}
  Renormalization group ideas and results from critical droplet theory are
  used to construct a scaling ansatz for the imaginary component of the free
  energy of an Ising model in its metastable state close to the critical
  point. The analytic properties of the free energy are used to determine
  scaling functions for the free energy in the vicinity of the critical point
  and the abrupt transition line. These functions have essential singularities
  at the abrupt transition. Analogous forms for the magnetization and susceptibility in
  two dimensions are fit to numeric data and show good agreement, especially
  when their nonsingular behavior is modified to match existing numeric results.
\end{abstract}

\maketitle

The Ising model is the canonical example of a system with a continuous phase
transition, and the study of its singular properties marked the first success
of the renormalization group (\textsc{rg}) method in statistical physics
\cite{wilson.1971.renormalization}. Its status makes sense: it's a simple
model whose continuous phase transition contains all the essential features of
more complex ones, but admits \textsc{rg} methods in a straightforward way and
has exact solutions in certain dimensions and for certain parameter
restrictions. However, the Ising critical point is not simply a continuous
transition: it also ends a line of abrupt phase transitions extending from it
at zero field below the critical temperature. Though typically neglected in
\textsc{rg} scaling analyses of the critical point, we demonstrate that there
are numerically measurable contributions to scaling due to the abrupt
transition line that cannot be accounted for by analytic changes of control or
thermodynamic variables.

\textsc{Rg} analysis predicts that the singular part of the free energy per
site $F$ as a function of reduced temperature $t=1-T_\c/T$ and field $h=H/T$
in the vicinity of the critical point takes the scaling form
$F(t,h)=|t|^{2-\alpha}\fF(h|t|^{-\beta\delta})$ for the low temperature phase
$t<0$ \cite{cardy.1996.scaling}. When studying the properties of the Ising
critical point, it is nearly always assumed that the universal scaling
function $\fF$ is analytic, i.e., has a convergent Taylor series. However, it
has long been known that there exists an essential singularity in $\fF$ at
zero argument, though its effects have long been believed to be unobservable
\cite{fisher.1967.condensation}, or simply just neglected
\cite{guida.1997.3dising,schofield.1969.parametric,schofield.1969.correlation,caselle.2001.critical,josephson.1969.equation,fisher.1999.trigonometric}.
With careful analysis, we have found that assuming the presence of the
essential singularity is predictive of the scaling form of, for instance, the
susceptibility and magnetization.

The provenance of the essential singularity can be understood using the
methods of critical droplet theory for the decay of an Ising system in a
metastable state, i.e., an equilibrium Ising state for $T<T_\c$, $H>0$
subjected to a small negative external field $H<0$. The existence of an
essential singularity has also been suggested by transfer matrix
\cite{mccraw.1978.metastability,enting.1980.investigation,mangazeev.2008.variational,mangazeev.2010.scaling} and \textsc{rg}
methods \cite{klein.1976.essential}, and a different kind of essential
singularity is known to exist in the zero-temperature susceptibility
\cite{orrick.2001.susceptibility,chan.2011.ising,guttmann.1996.solvability,nickel.1999.singularity,nickel.2000.addendum,assis.2017.analyticity}.  It has long been known that the decay
rate $\Gamma$ of metastable states in statistical mechanics is often related
to the metastable free energy $F$ by $\Gamma\propto\im F$
\cite{langer.1969.metastable,penrose.1987.rigorous,gaveau.1989.analytic,privman.1982.analytic}.
`Metastable free energy' can be thought of as either an analytic continuation
of the free energy through the abrupt phase transition, or restriction of the
partition function trace to states in the vicinity of the local free energy
minimum that characterizes the metastable state. In any case, the free energy
develops a nonzero imaginary part in the metastable region. Heuristically,
this can be thought of as similar to what happens in quantum mechanics with a
non-unitary Hamiltonian: the imaginary part describes loss rate of probability
that the system occupies any `accessible' state, which corresponds to decay.

In critical droplet theory, the metastable state decays when a domain of the
equilibrium state forms whose surface-energy cost for growth is outweighed by
bulk-energy gains. There is numerical evidence that, near the critical point,
these droplets are spherical \cite{gunther.1993.transfer-matrix}. The free energy
cost of the surface of a droplet of radius $R$ is $\Sigma S_\dim R^{\dim-1}$
and that of its bulk is $-M|H|V_\dim R^\dim$, where $S_\dim$ and $V_\dim$ are
the surface area and volume of a $(\dim-1)$-sphere, respectively, and $\Sigma$
is the surface tension of the equilibrium--metastable interface. The critical
droplet size is then $R_\c=(\dim-1)\Sigma/M|H|$ and the free energy of the
critical droplet is $\Delta F_\c=\pi^{\dim/2}\Sigma^\dim((\dim-1)/M|H|)^{\dim-1}/\Gamma(1+\dim/2)$.
Assuming the typical singular scaling forms
$\Sigma/T=|t|^\mu\fS(h|t|^{-\beta\delta})$ and $M=|t|^\beta\mathcal
M(h|t|^{-\beta\delta})$ and using known hyperscaling relations
\cite{widom.1981.interface}, this implies a scaling form
\def\eqcritformone{
  T\frac{\pi^{\dim/2}(\dim-1)^{\dim-1}}{\Gamma(1+\dim/2)}
    \frac{\fS^\dim(h|t|^{-\beta\delta})}{(-h|t|^{-\beta\delta}
    \fM(h|t|^{-\beta\delta}))^{\dim-1}}
}
\def\eqcritformtwo{
  T\fG^{-(\dim-1)}(h|t|^{-\beta\delta})
}
\ifreprint
\[
  \begin{aligned}
    \Delta F_\c
    &=\eqcritformone\\
    &\sim\eqcritformtwo
  \end{aligned}
\]
\else
\[
  \Delta F_\c=\eqcritformone\sim\eqcritformtwo
\]
\fi
for the free energy change due to the critical droplet.
Since both surface tension and magnetization are finite and nonzero for $H=0$
at $T<T_\c$, $\fG(X)=-BX+\O(X^2)$ for small negative $X$ with
\[
  B=\frac{\fM(0)}{\dim-1}\bigg(\frac{\Gamma(1+\dim/2)}
    {\pi^{\dim/2}\fS(0)^\dim}\bigg)^{1/(\dim-1)}.
\]
This first term in the scaling function $\fG$ is related to the ratio between
the correlation length $\xi$ and the critical domain radius $R_\c$, with
\[
  Bh|t|^{-\beta\delta}
  =\frac\xi{R_\c}\bigg(\frac{\Gamma(1+\dim/2)}
    {\pi^{\dim/2}\fS(0)(\xi_0^-)^{\dim-1}}\bigg)^{1/(\dim-1)}
\]
where the critical amplitude for the correlation length $\xi_0^-$ is defined
by $\xi=\xi_0^-|t|^{-\nu}$ for $t<T_\c$ and $H=0$. Since $\fS(0)(\xi_0^-)^{\dim-1}$ is a
universal amplitude ratio \cite{zinn.1996.universal},
$(Bh|t|^{-\beta\delta})/(\xi/R_\c)$ is a universal quantity.  The decay rate
of the metastable state is proportional to the Boltzmann factor corresponding
to the
creation of a critical droplet, yielding
\[
  \im F\sim\Gamma\propto e^{-\Delta F_\c/T}
    =e^{-\fG(h|t|^{-\beta\delta})^{-(\dim-1)}}.
\]
For $d>1$ this function has an essential singularity in the invariant
combination $h|t|^{-\beta\delta}$.

This form of $\im F$ for small $h$ is well known
\cite{langer.1967.condensation,harris.1984.metastability}.  We make the
scaling ansatz that the imaginary part of the metastable free energy has the
same singular behavior as the real part of the equilibrium free energy, and
that for small $t$, $h$, $\im F(t,h)=|t|^{2-\alpha}\fiF(h|t|^{-\beta\delta})$,
where
\[
  \fiF(X)=A\Theta(-X)(-BX)^be^{-1/(-BX)^{\dim-1}}
  \label{eq:im.scaling}
\]
and $\Theta$ is the Heaviside function. Results from combining an analysis of
fluctuations on the surface of critical droplets with \textsc{rg} recursion
relations suggest that $b=-(d-3)d/2$ for $d=2,4$ and $b=-7/3$ for $d=3$
\cite{houghton.1980.metastable,rudnick.1976.equations,gunther.1980.goldstone}.
Assuming that $F$ is analytic in the upper complex-$h$ plane, the real part of
$F$ in the equilibrium state can be extracted from this imaginary metastable
free energy using the Kramers--Kronig relation
\[
  \re F(t,h)=\frac1\pi\int_{-\infty}^\infty\frac{\im F(t,h')}{h'-h}\,\dd h'.
  \label{eq:kram-kron}
\]
This relationship has been used to compute high-order moments of the free
energy with $H$ in good agreement with transfer matrix expansions
\cite{lowe.1980.instantons}. Here, we evaluate the integral explicitly to come
to
functional forms.  In \threedee\ and \fourdee\ this can be done directly
given our scaling ansatz, yielding
\def\eqthreedeeone{
  \fF^\threedee(Y/B)&=
  \frac{A}{12}\frac{e^{-1/Y^2}}{Y^2}
  \bigg[\Gamma(\tfrac16)E_{7/6}(-Y^{-2})
}
\def\eqthreedeetwo{
  -4Y\Gamma(\tfrac23)E_{5/3}(-Y^{-2})\bigg]
}
\def\eqfourdeeone{
  \fF^\fourdee(Y/B)&=
  -\frac{A}{9\pi}\frac{e^{1/Y^3}}{Y^2}
  \Big[3\ei(-Y^{-3})
}
\def\eqfourdeetwo{
  +3\Gamma(\tfrac23)\Gamma(\tfrac13,Y^{-3})
  +\Gamma(\tfrac13)\Gamma(-\tfrac13,Y^{-3})\Big]
}
\ifreprint
\begin{align}
  &\begin{aligned}
    \eqthreedeeone\\
    &\hspace{7em}
    \eqthreedeetwo
  \end{aligned}
  \\
  &\begin{aligned}
    \eqfourdeeone
    \\
    &\hspace{-0.5em}
    \eqfourdeetwo,
  \end{aligned}
\end{align}
\else
\begin{align}
  \eqthreedeeone\eqthreedeetwo
  \\
  \eqfourdeeone\eqfourdeetwo,
\end{align}
\fi
where $E_n$ is the generalized exponential integral and $\Gamma(x,y)$ is the
incomplete gamma function.
At the level of truncation of \eqref{eq:im.scaling} at which we are working
the Kramers--Kronig relation does not converge in \twodee. However, higher
moments can still be extracted, e.g., the susceptibility, by taking
\[
  \chi=\pd MH=-\frac1{T}\pd[2]Fh
  =-\frac2{\pi T}\int_{-\infty}^\infty\frac{\im F(t,h')}{(h'-h)^3}\,\dd h'.
\]
With a scaling form defined by $T\chi=|t|^{-\gamma}\fX(h|t|^{-\beta\delta})$,
this yields
\[
  \fX^\twodee(Y/B)=\frac{AB^2}{\pi Y^3}\big[Y(Y-1)-e^{1/Y}\ei(-1/Y)\big]
  \label{eq:sus_scaling}
\]
Scaling forms for the free energy can then be extracted by direct integration
and their constants of integration fixed by known zero-field values, yielding
\begin{align}
  \label{eq:mag_scaling}
  \fM^\twodee(Y/B)
    &=\fM(0)+\frac{AB}{\pi}\bigg(1-\frac{Y-1}Ye^{1/Y}\ei(-1/Y)\bigg)\\
  \fF^\twodee(Y/B)
    &=-Y\bigg(\frac{\fM(0)}B-\frac{A}\pi e^{1/Y}\ei(-1/Y)\bigg)
    \label{eq:2d_free_scaling}
\end{align}
with $F(t,h)=t^2\fF(h|t|^{-15/8})+t^2\log t^2$ in
two dimensions, as $\alpha=0$ and $\beta\delta=\frac{15}8$.

How are these functional forms to be interpreted? Though the scaling function
\eqref{eq:im.scaling} for the imaginary free energy of the metastable state is
asymptotically correct sufficiently close to the critical point, the results
of the integral relation \eqref{eq:kram-kron} are not, since there is no limit
of $t$ or $h$ in which it becomes arbitrarily correct for a given truncation
of \eqref{eq:im.scaling}. It is well established that this method of using
unphysical or metastable elements of a theory to extract properties of the
stable or equilibrium theory is only accurate for high moments of those
predictions \cite{parisi.1977.asymptotic,bogomolny.1977.dispersion}.  The
functions above should be understood as possessing exactly the correct
singularity  at the coexistence line, but requiring polynomial corrections,
especially for smaller integer powers. Using these forms in conjunction with
existing methods of describing the critical equation of state or critical
properties with analytic functions in $h$ will incorporate these low-order
corrections while preserving the correct singular structure. In other words,
the scaling functions can be \emph{exactly} described by
$\tilde\fF(X)=\fF(X)+f(X)$ for some analytic function $f$. Higher order terms
in the expansion of $\tilde\fF$ become asymptotically equal to those of $\fF$
because, as an analytic function, progressively higher order terms of $f$ must
eventually become arbitrarily small \cite{flanigan.1972.complex}.

How predictive are these scaling forms in the proximity of the critical point
and the abrupt transition line? We simulated the \twodee\ Ising model on square
lattice using a form of the Wolff algorithm modified to remain efficient in
the presence of an external field. Briefly, the external field $H$ is applied
by adding an extra spin $s_0$ with coupling $|H|$ to all others
\cite{dimitrovic.1991.finite}. A quickly converging estimate for the
magnetization in the finite-size system was then made by taking
$M=\sgn(H)s_0\sum s_i$, i.e., the magnetization relative to the external spin
\cite{kent-dobias.2018.wolff}.
Data was then taken for susceptibility and magnetization for
$T_\c-T,H\leq0.1$. This data, rescaled as appropriate to collapse onto a
single curve, is plotted in Fig.~\ref{fig:scaling_fits}.

For the \twodee\ Ising model on a square lattice, exact results at zero
temperature have $\fS(0)=4/T_\c$, $\fM(0)=(2^{5/2}\arcsinh1)^\beta$
\cite{onsager.1944.crystal}, and $\fX(0)=C_0^-=0.025\,536\,971\,9$ \cite{barouch.1973.susceptibility}, so that
$B=T_\c^2\fM(0)/\pi\fS(0)^2=(2^{27/16}\pi(\arcsinh1)^{15/8})^{-1}$. If we
assume incorrectly that \eqref{eq:sus_scaling} is the true asymptotic form of
the susceptibility scaling function, then
$T\chi(t,0)|t|^\gamma=\lim_{X\to0}\fX^\twodee(X)=2AB^2/\pi$ and the constant
$A$ is fixed to $A=\pi\fX(0)/2B^2=2^{19/8}\pi^3(\arcsinh1)^{15/4}C_0^-$.  The
resulting scaling functions $\fX$ and $\fM$ are plotted as solid blue lines in
Fig.~\ref{fig:scaling_fits}. Though there is good agreement
between our functional forms and what is measured, there
are systematic differences that can be seen most clearly in the
magnetization. This is to be expected based on our earlier discussion: these
scaling forms should only be expected to well-describe the singularity at the
abrupt transition. Our forms both exhibit incorrect low-order
coefficients at the transition (Fig.~\ref{fig:series}) and incorrect
asymptotics as $h|t|^{-\beta\delta}$ becomes very large.

In forthcoming work,
we develop a method to incorporate the essential singularity in the scaling
functions into a form that also incorporates known properties of the scaling
functions in the rest of the configuration space using a Schofield-like
parameterization \cite{schofield.1969.parametric,caselle.2001.critical,kent-dobias.2018.parametric}. Briefly, we define
parameters $R$ and $\theta$ by
\begin{align}
  t=R(1-\theta^2)
  &&
  h=h_0R^{\beta\delta}g(\theta)
\end{align}
where $h_0$ is a constant and $f$ is an arbitrary odd function whose first
finite zero $\theta_\c>1$ corresponds to the abrupt transition. In these
coordinates the invariant combination $h|t|^{-\beta\delta}$ is given by
\[
  h|t|^{-\beta\delta}=\frac{h_0g(\theta)}{|1-\theta^2|^{\beta\delta}}=\frac{h_0(-g'(\theta_\c))}{(\theta_c^2-1)^{\beta\delta}}(\theta_\c-\theta)
  +\O\big( (\theta_\c-\theta)^2\big),
\]
an analytic function of $\theta$ about $\theta_\c$.
The simplest
function of the coordinate $\theta$ that exhibits the correct singularity at
the abrupt transitions at $h=\pm0$, $t<0$ is
\[
  \fX(\theta)=\fX\bigg(\frac{h_0(-g'(\theta_\c))}{(\theta_\c^2-1)^{\beta\delta}}(\theta_\c-\theta)\bigg)+\fX\bigg(\frac{h_0(-g'(\theta_\c))}{(\theta_\c^2-1)^{\beta\delta}}(\theta_\c+\theta)\bigg)
\]
This function is analytic in the range $-\theta_\c<\theta<\theta_\c$.
In order to correct its low-order behavior to match that expected, we
make use of both the freedom of the coordinate transformation $g$ and an
arbitrary analytic additive function $Y$,
\begin{align}
  g(\theta)=\bigg(1-\frac{\theta^2}{\theta_\c^2}\bigg)\sum_{n=0}^\infty
  g_n\theta^{2n+1}
  &&
  Y(\theta)=\sum_{n=0}^\infty Y_n\theta^{2n}
\end{align}
so that $\tilde\fX(\theta)=\fX(\theta)+Y(\theta)$. By manipulating these
coefficients, we can attempt to give the resulting scaling form a series
expansion consistent with known values. One such prediction---made by fixing
the first four terms in the low-temperature, critical isotherm, and
high-temperature expansions 
of $\tilde\fX$---is shown as a dashed yellow line in
Fig.~\ref{fig:scaling_fits}.
As shown in Fig.~\ref{fig:series}, the low-order free energy coefficients of this prediction match known values
exactly up to $n=5$, and improve the agreement with higher-order coefficients.
Unlike scaling forms which treat $\fX$ as analytic at the coexistence line,
the series coefficients of the scaling form developed here increase without
bound at high order.


\begin{figure}
  \input{fig-susmag}
  \caption{
    Scaling functions for (top) the susceptibility and (bottom) the
    magnetization plotted in terms of the invariant combination
    $h|t|^{-\beta\delta}$. Points with error bars show data with sampling
    error taken from simulations of a $4096\times4096$ square-lattice Ising
    model with periodic boundary conditions and $T_\c-T=0.01,0.02,\ldots,0.1$
    and $H=0.1\times(1,2^{-1/4},\ldots,2^{-50/4})$. The solid blue lines
    show our analytic results \eqref{eq:sus_scaling} and
    \eqref{eq:mag_scaling}, the dashed yellow lines show 
    a scaling function modified to match known series expansions of the
    susceptibility
    to third order, and the
    dotted green lines show the
    polynomial resulting from truncating the known series expansion after the eight terms
    reported by \cite{mangazeev.2008.variational,mangazeev.2010.scaling}.
  }
  \label{fig:scaling_fits}
\end{figure}

\begin{figure}
  \input{fig-series}
  \caption{
    The series coefficients defined by $\tilde\fF(X)=\sum_nf_nX^n$. The blue
    pluses correspond to the scaling form \eqref{eq:2d_free_scaling}, the
    yellow saltires correspond to a scaling function modified to match known
    series expansions of the susceptibility to third order---and therefore
    the free energy to fifth order---and the green
    stars
    correspond to the first eight coefficients from
    \cite{mangazeev.2008.variational,mangazeev.2010.scaling}. The modified
    scaling function and the known coefficients match exactly up to $n=5$.
  }
  \label{fig:series}
\end{figure}

Abrupt phase transitions, such as the jump in magnetization in the Ising
model below $T_\c$, are known to imply essential singularities in the free
energy that are usually thought to be unobservable in practice. We have
argued that this essential singularity controls the universal scaling
behavior near continuous phase transitions, and have derived an explicit
analytical form for the singularity in the free energy, magnetization,
and susceptibility for the Ising model. We have developed a Wolff algorithm 
for the Ising model in a field, and showed that incorporating our singularity
into the scaling function gives good convergence to the simulations in \twodee.

Our results should allow improved high-precision functional forms for the free
energy~\cite{caselle.2001.critical}, and should have implications for the scaling
of correlation functions~\cite{chen.2013.universal,wu.1976.spin}. Our methods might be generalized
to predict similar singularities in systems where nucleation and metastability
are proximate to continuous phase transitions, such as 2D superfluid
transitions~\cite{ambegaokar.1978.dissipation,ambegaokar.1980.dynamics}, the melting of 2D crystals~\cite{dahm.1989.dynamics}, and
freezing transitions in glasses, spin glasses, and other disordered systems.


%We have used results from the properties of the metastable Ising ferromagnet
%and the analytic nature of the free energy to derive universal scaling
%functions for the free energy, and in \twodee the magnetization and
%susceptibility, in the limit of small $t<0$ and $h$. Because of an essential
%singularity in these functions at $h=0$---the abrupt transition line---their
%form cannot be brought into that of any regular function by analytic
%redefinition of control or thermodynamic variables. These predictions match
%the results of simulations well. Having demonstrated that the essential
%singularity in thermodynamic functions at the abrupt transition leads to
%observable scaling effects, we hope that these functional forms will be used in
%conjunction with traditional perturbation methods to better express the
%equation of state of the Ising model in the whole of its parameter space.

\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{essential-ising}

\end{document}