# Slow quench dynamics in classical systems: kinetic Ising model and zero-
range process

Priyanka1, Sayani Chatterjee2,†, Kavita Jain2 1 Department of Physics (MC
0435) and Center for Soft Matter and Biological Physics, Virginia Tech,
Robeson Hall, 850 West Campus Drive, Blacksburg, VA 24061, USA
2Theoretical Sciences Unit, Jawaharlal Nehru Centre for Advanced Scientific
Research, Bangalore 560064, India
† Present address: No. 19 Railway Gate, P. O. Bengal Enamel, North 24
Parganas, West Bengal 743122, India

###### Abstract

While a large number of studies have focused on the nonequilibrium dynamics of
a system when it is quenched instantaneously from a disordered phase to an
ordered phase, such dynamics have been relatively less explored when the
quench occurs at a finite rate. Here we study the slow quench dynamics in two
paradigmatic models of classical phase transitions, viz., one-dimensional
kinetic Ising model and mean-field zero-range process, when the system is
annealed slowly to the critical point. Starting from the time evolution
equations for the spin-spin correlation function in the Ising model and mass
distribution in the zero-range process, we derive the Kibble-Zurek scaling
laws. We then test a recent proposal that critical coarsening which is ignored
in the Kibble-Zurek argument plays a role in the nonequilibrium dynamics close
to the critical point. We find that the defect density in the Ising model and
scaled mass distribution in the zero-range process decay linearly with the
time remaining until the end of quench provided the final quench point is
approached sufficiently fast, and sublinearly otherwise. As the linear scaling
also holds when a system following an instantaneous quench is allowed to
coarsen for a finite time interval, we conclude that critical coarsening
describes the slow quench dynamics in the vicinity of the critical point only
if the annealing is not too slow.

Keywords: coarsening processes, kinetic Ising models, zero-range processes

##  1 Introduction

The nonequilibrium dynamics following rapid quenches where the annealing time
from a disordered phase to an ordered phase is much shorter than the time
scale over which phase ordering occurs have been extensively studied [1]. In
the recent past, attention has focused on slow quenches that are relevant to
understand the residual density of defect structures in the slowly cooling
early universe [2], glassy states that are obtained on cooling a liquid at a
finite rate [3, 4] and more generally, the nonequilibrium dynamics of
classical [5, 6, 7, 8, 9, 10, 11, 12, 13] and quantum systems [14, 15, 16, 17,
18, 19] that exhibit second order phase transitions [20]. Theoretical
predictions for quantum annealed systems have also been tested in experiments
with different media and especially, ultra-cold atomic gases [21, 22]. A
general conclusion of these investigations is that there are more defects at
the end of slow quench than in equilibrium which may be understood using an
argument based on equilibrium state properties [2, 20, 23]. The idea is that
as the correlation length is of order unity in the disordered phase, the
system can equilibrate to the changing temperature (or the relevant tuning
parameter). However, close to the critical point where the correlation length
diverges, the system falls out of equilibrium and evolves very slowly. Then,
assuming that the system does not evolve at all in the nonequilibrium regime,
the scaling of the density of excess defects at the end of quench with
annealing rate can be predicted and has been verified in several studies [24,
25].

The ‘frozen dynamics’ assumption has been questioned recently, and using
scaling ideas and numerical simulations, it has been argued that the
nonequilibrium regime is characterized by critical coarsening which can lead
to scaling laws different from the Kibble-Zurek prediction [6]. To the best of
our knowledge, this proposal has not been investigated analytically;
furthermore, most studies have focused on the properties at the end of the
quench and are limited to the linear cooling of the system. Here, we study the
one-dimensional Ising model with Glauber dynamics for general cooling schemes
analytically, and derive scaling properties of the defect density at and close
to the critical point. We then consider a zero-range process with time-
dependent hop rates and study these properties both analytically and
numerically.

It is well known that the classical ferromagnetic Ising model exhibits a
nontrivial phase transition above one dimension. The phase ordering dynamics
on slowly cooling the Ising system from a disordered phase to the critical
point or ordered phase have been studied, mainly numerically, in $d\geq 2$
dimensions without [6, 8, 26, 10, 11] and with [27, 28, 29, 30] quenched
disorder. In Sec. 2, we study the one-dimensional pure Ising model with
Glauber dynamics as it is analytically tractable. Using the known exact
solution for the spin-spin correlation function with time-dependent
temperature [31, 32, 33], we calculate the defect density at the end of the
quench. Our expression for the dependence of residual defect density on
cooling rate matches with the predictions of the Kibble-Zurek argument [2, 20,
23] as well as the previous work on one-dimensional Ising model [33, 7]. We
also calculate the prefactor for the excess defect density (note that the
prefactors calculated in [7] are incorrect, as shown here), and find that the
density of defects is substantially smaller than that predicted under ‘frozen
dynamics’ assumption. We then leverage the exact solution to investigate the
role of critical coarsening in the decrease of the defect density [6], and
find that for sufficiently fast cooling (at finite rate), the defect density
decays with the remaining time until the end of quench linearly and
sublinearly otherwise. As the scaling in the former case is the same as that
when the system is quenched instantaneously and allowed to coarsen for a
finite time, our main conclusion is that the dynamics in the nonequilibrium
regime can be described by the instantaneous quench model, provided the final
temperature is approached fast enough.

In Sec. 3, we consider a zero-range process in mean field geometry and exploit
the insights gained from the Ising model to understand its slow quench
dynamics. The stationary state of this model is known exactly in arbitrary
dimensions and for hop rates considered here, it exhibits a phase transition
from a fluid phase with an order unity particles distributed per site to a
condensate phase where a macroscopic number of particles occupy a site, as the
parameter $b$ in the hop rate is increased [34]. The instantaneous quench
dynamics have been studied using scaling arguments in mean-field geometry and
in one dimension [35], and the slow quench dynamics in the one-dimensional
zero-range process were studied numerically in [9]. Here, we develop a scaling
theory for the zero-range process with time-dependent rates when it is
annealed slowly to the critical point and derive the Kibble-Zurek scaling laws
for the mass distribution. Our differential equation for the scaling function
does not appear to be solvable, and we therefore study the mass distribution
close to the critical point numerically. Analogous to the Ising model, we find
that it varies linearly with the remaining time, provided the parameter $b$
approaches the critical point sufficiently fast.

In Sec. 4, we summarize our results and some open questions.

##  2 One-dimensional Ising model with Glauber dynamics

We consider a one-dimensional Ising model on a ring with $L$ sites. In the
absence of external field, the Hamiltonian
$H=-J\sum_{i=1}^{L}\sigma_{i}\sigma_{i+1}$ where the Ising spin
$\sigma_{i}=\pm 1$ and the coupling $J>0$. Due to periodic boundary condition,
the spin variable $\sigma_{L+1}=\sigma_{1}$. Under Glauber dynamics, a spin
configuration $\\{\sigma_{1},...,\sigma_{j},...\sigma_{L}\\}$ evolves via
single-spin flips, and the distribution of the configuration obeys the
following master equation [36]

$\frac{\partial\psi(...,\sigma_{j},...)}{\partial
t}=\sum_{j}W(-\sigma_{j}\to\sigma_{j}){\psi}(...,-\sigma_{j},...)-W(\sigma_{j}\to-\sigma_{j}){\psi}(...,\sigma_{j},...)~{},$
(1)

where the transition rate

$W(\sigma_{j}\to-\sigma_{j})=\frac{A}{2}\left(1-\frac{\gamma}{2}\sigma_{j}(\sigma_{j-1}+\sigma_{j+1})\right)~{},$
(2)

and

$\displaystyle\gamma(T)$ $\displaystyle=$
$\displaystyle\tanh\left(\frac{2J}{k_{B}T}\right)~{},$ (3) $\displaystyle
A(T)$ $\displaystyle=$ $\displaystyle
A_{0}\exp\left(-\frac{\Delta}{k_{B}T}\right)~{}.$ (4)

As the parameter $\gamma$ decreases with temperature $T$, by virtue of (2),
the spins have a tendency to align with each other at low temperatures. The
Arrhenius rate $A$ introduces an activation energy $\Delta$ for a spin flip to
occur and is evidently important at low temperatures [31, 33].

Motivated by the glass-problem in which the dynamics are ‘frozen’ in a slowly
cooled liquid at low temperatures due to the activation barrier [4], the slow
quench dynamics in the one-dimensional Ising model have been studied
analytically by including a time-dependence in the temperature [31, 32, 28,
33]. However, even in the absence of an activation barrier, the dynamics slow
down in the vicinity of the critical point as the relaxation time becomes
longer than the annealing time [2, 20, 23]; here, we are interested in
understanding such dynamics and therefore, we will set $\Delta=0$ throughout
the following discussion. For convenience, we will also set $J/k_{B}=1$.

###  2.1 General expression for two-point correlation function

We consider the two-point correlation function
$G_{k}(t)=\langle\sigma_{i}\sigma_{i+k}\rangle$ where the angular brackets
denote an average with respect to the ensemble distribution
${\psi}(\vec{\sigma})$ (for a generalization to higher-point correlation
function, see [37]). For arbitrary $\gamma$, it obeys the following exact
equation [36, 31],

$\displaystyle{\dot{G}_{k}}=-2G_{k}+\gamma(t)(G_{k-1}+G_{k+1})~{},~{}k=1,...,L-1~{},$
(5)

where the dot denotes a time derivative and the boundary condition
$G_{0}(t)=G_{L}(t)=1$.

When the temperature is constant in time, the dynamics of the two-point
correlation function have been studied in detail [36] (see Cornell in [38] for
a review). On an infinite ring, the equilibrium correlation function is given
by

$G_{k}(t\to\infty)\equiv{\cal
G}_{k}=\left(\frac{1-\sqrt{1-\gamma^{2}}}{\gamma}\right)^{k}~{},~{}k=0,1,...$
(6)

where we have used that the correlation function is bounded above by one. In
the stationary state, close to the critical point, the relaxation time scales
with the correlation length with the dynamic critical exponent $z=2$. If the
system is instantaneously quenched from high temperature to a low temperature,
the equal time correlation function ${\hat{\cal G}}_{k}(t)={\hat{\cal
F}}(k/\sqrt{t})$ which thus yields the coarsening exponent ${\hat{z}}=2$.

We are interested in understanding the dynamics when the Ising system
initially in an equilibrium state at infinite temperature is cooled to zero
temperature in finite time $\tau$ (the assumption of infinite $T(0)$ is not
especially restrictive but simplifies the analysis). For time-dependent
$\gamma$, an explicit expression for the two-point correlation function has
been obtained in [31] and [33]. Here, we follow the latter where it has been
shown that on an infinite ring, the deviation of the correlation function
$G_{k}(t)$ at time $t$ from the corresponding equilibrium value at
instantaneous temperature is given by [33]

$\displaystyle{\cal G}_{k}(t)-G_{k}(t)$ $\displaystyle=$
$\displaystyle\int_{0}^{\pi}dq\phi_{k}(q)\int_{0}^{t}dt^{\prime}e^{-2\int_{t^{\prime}}^{t}dy(1-\gamma(y)\cos
q)}\sum_{m=1}^{\infty}\phi_{m}(q){\dot{\cal G}}_{m}(t^{\prime})$ (7)
$\displaystyle+$
$\displaystyle\int_{0}^{\pi}dq\phi_{k}(q)e^{-2\int_{0}^{t}dy(1-\gamma(y)\cos
q)}\sum_{m=1}^{\infty}\phi_{m}(q)({\cal G}_{m}(0)-G_{m}(0))~{},$

where $\phi_{k}(q)=\sqrt{\frac{2}{\pi}}\sin(kq)$ obeys the eigenvalue equation
$\gamma(\phi_{k-1}(q)+\phi_{k-1}(q))-2\phi_{k}(q)=-\lambda\phi_{k}(q)$ with
eigenvalue $\lambda(q)=2(1-\gamma\cos q)$. The above expression is obtained on
expanding ${\cal G}_{k}(t)-G_{k}(t)$ as a linear combination of the
eigenfunctions $\phi_{k}(q)$ and using the time evolution equation (5).
Furthermore, since the system is assumed to be initially in the equilibrium
state at infinite temperature so that $G_{k}(0)={\cal G}_{k}(0)=\delta_{k,0}$,
on carrying out the sum over $m$ on the right-hand side (RHS) of (7) exactly,
we obtain

$\displaystyle{\cal
G}_{k}(t)-G_{k}(t)=\int_{0}^{t}dt^{\prime}e^{-2(t-t^{\prime})}{\dot{\gamma}(t^{\prime})}\int_{0}^{\pi}\frac{dq}{\pi}\frac{\sin(q)\sin(kq)}{(1-\gamma(t^{\prime})\cos(q))^{2}}e^{2\cos
q\int_{t^{\prime}}^{t}dy\gamma(y)}~{}.$ (8)

In the following subsections, we will analyze the above double integral in
detail.

Following previous work [33, 7], we consider the following class of cooling
protocols,

$\displaystyle\gamma(t,\tau)=$ $\displaystyle 1-\left(1-x\right)^{\alpha}$ ,
$\alpha>0$ (9a) $\displaystyle\gamma(t,\tau)=$ $\displaystyle
1-\exp\left[1-\left(1-x\right)^{-\beta}\right]$ , $\beta>0$ , (9b)

where $x=t/\tau$. The parameter $\gamma(x)$ increases monotonically from zero
to one with increasing time and, (9a) and (9b), respectively, correspond to
temperature,

$\displaystyle T(x)\stackrel{{\scriptstyle t\to\tau}}{{\sim}}$
$\displaystyle\frac{-4}{\alpha\ln(1-x)}$ (logarithmic cooling) (10a)
$\displaystyle T(x)\stackrel{{\scriptstyle t\to\tau}}{{\sim}}$ $\displaystyle
4\left(1-x\right)^{\beta}$ (algebraic cooling) . (10b)

Note that the above cooling schemes reduce to the instantaneous quench limit
when $\alpha,\beta\to\infty$.

###  2.2 Dynamics deep in the paramagnetic phase

We first consider the parameter regime where $\gamma\ll 1$. Since the
integrand in the outer integral of (8) suppresses large $t-t^{\prime}$, the
inner integral receives most contribution when $\gamma$ is evaluated in the
vicinity of $t$ thus yielding

$\displaystyle{\cal G}_{k}(t)-G_{k}(t)$ $\displaystyle\approx$
$\displaystyle\int_{0}^{t}dt^{\prime}\int_{0}^{\pi}\frac{dq}{\pi}\frac{\dot{\gamma}(t)e^{-2(t-t^{\prime})(1-\gamma(t)\cos(q))}\sin(q)\sin(kq)}{(1-\gamma(t)\cos(q))^{2}}$
(11) $\displaystyle=$
$\displaystyle\frac{{\dot{\gamma}(t)}}{2\pi}\int_{0}^{\pi}dq\frac{\sin(kq)\sin(q)}{(1-\gamma(t)\cos(q))^{3}}$
(12) $\displaystyle=$
$\displaystyle\frac{k}{4}\frac{{\dot{\gamma}}}{\gamma}\frac{{\cal
G}_{k}(t)}{1-\gamma^{2}}\left(k+\frac{1}{\sqrt{1-\gamma^{2}}}\right)~{}.$ (13)

At short times, as Fig. 1 shows, the correlation function ${G}_{k}(t)$ is
close to the equilibrium correlation function ${\cal G}_{k}(T(t))$ at
instantaneous temperature. But, with increasing time, the deviation increases
as the system relaxes slower than the rate of change in temperature. For a
given $t\ll\tau$, the deviation ${\cal G}_{k}(t)-G_{k}(t)\propto k^{2}x^{k}$
is, however, nonmonotonic in the distance between the spins with the peak
occurring at larger $k$ for larger times. Figure 1 also shows that the
nonequilibrium correlation function, $G_{k}(t)$ is always smaller than the
equilibrium correlation function, ${\cal G}_{k}(t)$; this is because the
nonequilibrium system has not completely relaxed to the stationary state at
the instantaneous temperature and is effectively at a slightly higher
temperature.

###  2.3 Dynamics close to the critical point

Close to zero temperature, the behavior of the correlation function is
determined by how the function $\gamma$ approaches one. We therefore consider
the logarithmic and algebraic cooling schemes here separately.

####  2.3.1 Logarithmic cooling

Using the expression (9a) for $\gamma$ in (8), we obtain

${\cal
G}_{k}(t)-G_{k}(t)=\frac{\alpha}{\pi}\int_{0}^{x}dy(1-x+y)^{\alpha-1}\int_{0}^{\pi}dq\frac{\sin(q)\sin(kq)e^{-2\tau
z(x,y,q)}}{\left[1-\cos(q)+(1-x+y)^{\alpha}\cos(q)\right]^{2}}~{},$ (14)

where

$z=(1-\cos(q))y+\left[\frac{(1-x+y)^{\alpha+1}-(1-x)^{\alpha+1}}{\alpha+1}\right]\cos(q)~{},$
(15)

and $x=t/\tau$. For $\tau\gg 1$ and $x\to 1$, the exponential term in the
inner integral on the RHS of (14) contributes when $q$ and $y$ are close to
zero. Therefore, for large $\tau$, we obtain

${\cal
G}_{k}(t)-G_{k}(t)\approx\frac{\alpha}{\pi}\int_{0}^{x}dy(1-x+y)^{\alpha-1}\int_{0}^{\infty}dq\frac{q\sin(kq)e^{-2\tau\left[\frac{yq^{2}}{2}+\frac{(1-x+y)^{\alpha+1}-(1-x)^{\alpha+1}}{\alpha+1}\right]}}{\left[\frac{q^{2}}{2}+(1-x+y)^{\alpha}\right]^{2}}~{}.$
(16)

Then a change of variables to
$Y=y\tau^{\frac{1}{1+\alpha}},Q=q\tau^{\frac{\alpha}{2(1+\alpha)}}$
immediately shows that the deviation in the correlation function has a scaling
form where the scaling function is given by

$\displaystyle{\cal G}_{k}(t)-G_{k}(t)\equiv F(K,Z)$ (17) $\displaystyle=$
$\displaystyle\frac{\alpha}{\pi}\int_{0}^{\infty}dY(Z+Y)^{\alpha-1}e^{-\frac{2}{\alpha+1}\left[(Z+Y)^{\alpha+1}-Z^{\alpha+1}\right]}\int_{0}^{\infty}dQ\frac{Q\sin(KQ)e^{-YQ^{2}}}{\left[\frac{Q^{2}}{2}+(Z+Y)^{\alpha}\right]^{2}}$
$\displaystyle=$
$\displaystyle\frac{1}{\pi}\int_{0}^{\infty}dY\exp\left[-\frac{2Y^{\frac{\alpha+1}{\alpha}}}{(\alpha+1)K^{\frac{2(\alpha+1)}{\alpha}}}\right]\int_{0}^{\infty}dQ^{\prime}\frac{Q^{\prime}\sin(Q^{\prime})e^{-\frac{Y^{1/\alpha}Q^{\prime
2}}{K^{2+(2/\alpha)}}}}{\left[\frac{Q^{\prime 2}}{2}+Y\right]^{2}}~{},$ (18)

and the scaling variables are

$K=k\tau^{-\frac{\alpha}{2(1+\alpha)}},Z=\frac{1-x}{1-{x_{*}}}=(1-x)\tau^{\frac{1}{1+\alpha}}~{}.$
(19)

Close to the critical point, since the equilibrium correlation function given
by (6) can be written as

${\cal G}_{k}(t)\approx e^{-k\sqrt{2(1-\gamma)}}=e^{-\sqrt{2}KZ^{\alpha/2}},$
(20)

the correlation function $G_{k}(t)$ is also a function of $K$ and $Z$. The
data collapse shown in Fig. 2 confirms that the correlation function
$G_{k}(\tau)$ is of the above scaling form at the critical point.

The scaling exponents in (19) are in agreement with the Kibble-Zurek argument
[2, 20, 23] which states that for $x<{x_{*}}$, the system stays close to the
equilibrium state at the instantaneous temperature. But, for $x>{x_{*}}$, it
falls out of equilibrium as the relaxation time exceeds the time remaining to
reach the final quench time. Therefore, at $x=x^{*}$, one obtains
$\tau(1-{x_{*}})\sim\xi^{z}$ where the correlation length
$\xi=2^{-1/2}(1-x)^{-\alpha/2}$ (see (20)) and the dynamic exponent $z=2$
[36]. This immediately yields the relevant time and length scale,
$1-{x_{*}}\sim\tau^{-\frac{1}{1+\alpha}}$ and
$\xi\sim\tau^{\frac{\alpha}{2(1+\alpha)}}$, respectively.

The double integral in (17) does not appear to be exactly solvable. But one
can find the scaling function $F(K,Z)$ in various limits of interest. We first
consider the deviation in the correlation function at the end of quench
($Z=0$), and find that

$\displaystyle 1-G_{k}(\tau)=F(K,0)=$
$\displaystyle\frac{K}{\sqrt{\pi}}\left(\frac{2}{1+\alpha}\right)^{\frac{1}{2(1+\alpha)}}\Gamma\left(\frac{1+2\alpha}{2+2\alpha}\right)$
, $K\ll 1$ (21a) $\displaystyle 1-G_{k}(\tau)=F(K,0)=$ $\displaystyle 1$ ,
$K\gg 1$ . (21b)

For finite $k$, on expanding the sine function in (16) for small argument and
performing the resulting integrals, we obtain (21a), and on taking the limit
$K\to\infty$ in (18), we arrive at (21b). The above expressions match with the
results obtained using numerical integration of (5) as shown in Fig. 2.

The expression (21a) can be used to find the domain wall density,
${D}(\tau)=(1-G_{1}(\tau))/2$ at the end of the quench; on comparing our
result for $D(\tau)$ with that in [7], we find that although the expression
(13) in [7] captures the correct $\tau$-dependence for the decay of domain
wall density, it underestimates the prefactor. Figure 3 shows that our (21a)
is in excellent agreement with the numerical data except for small $\alpha$;
this is because the subleading correction to (21a) is of order
$\tau^{-\alpha/(1+\alpha)}$, as can be seen by retaining the next order term
in (15) and using (19). Equation (21a) shows that the defect density at the
critical point decreases with increasing $\tau$, as intuitively expected. The
defect density $D(\tau)$ also decreases with increasing $\alpha$ which can be
understood if we reject the ‘frozen dynamics’ assumption [24, 25]: for large
$\alpha$, the system reaches low temperature quickly and spends rest of the
time equilibrating thereby decreasing the number of defects, while the reverse
holds for small $\alpha$. We also verify that in the instantaneous quench
limit $\alpha\to\infty$, (21a) reduces to (30) at $Z=0$ (see below).

####  2.3.2 Algebraic cooling

Here, we focus on the small-$k$ behavior. Using the cooling scheme (9b) in
(8), we find that

$\displaystyle{\cal G}_{k}(t)-G_{k}(t)$ (22) $\displaystyle\approx$
$\displaystyle\frac{k\beta}{\pi}\int_{0}^{x}dye^{1-(1-x+y)^{-\beta}}(1-x+y)^{-1-\beta}\int_{0}^{\pi}dq\frac{q\sin(q)e^{-2\tau
z(x,y,q)}}{[1+(e^{1-(1-x+y)^{-\beta}}-1)\cos(q)]^{2}}~{},$ (23)

where
$z=y(1-\cos(q))+\cos(q)\int_{x-y}^{x}dy^{\prime}e^{1-(1-y^{\prime})^{-\beta}}$.
As in the last subsection, we expand the integrand in the inner integral for
small $q$. Then, for large $\tau$ and $x\to 1$, we obtain

$\displaystyle{\cal G}_{k}(t)-G_{k}(t)$ $\displaystyle\approx$
$\displaystyle\frac{k}{\pi\sqrt{\tau}}\int_{0}^{\infty}\frac{dW}{\sqrt{y(W)}}\int_{0}^{\infty}dQ\frac{Q^{2}e^{-Q^{2}}}{(\frac{Q^{2}}{2}+W)^{2}}~{},$
(24)

where $W=\tau ye^{1-(1-x+y)^{-\beta}}$ is finite when $x\to 1,y\to
0,\tau\to\infty$. This implies that $1-x+y\approx(\ln\tau)^{-1/\beta}$ using
which we find that

$\displaystyle{\cal G}_{k}(t)-G_{k}(t)=$
$\displaystyle\frac{K}{\sqrt{\pi(1-Z)}}~{}$ , $Z<1$ (25a) $\displaystyle{\cal
G}_{k}(t)-G_{k}(t)=$ $\displaystyle 0~{}$ , $Z>1$ (25b)

where

$K=k\sqrt{\frac{(\ln\tau)^{1/\beta}}{\tau}}~{},~{}Z=(1-x)(\ln\tau)^{1/\beta}~{}.$
(26)

Equation (25a) shows that the deviation in the correlation function when
expressed in terms of scaling variables is independent of quench exponent
$\beta$, and equals (30) at $Z=0$ which, as mentioned above, is obtained when
$\alpha\to\infty$. This is a simple consequence of the fact that for $x\to 1$,
the function $\gamma(x)$ in (9b) decays faster than any power law. It can be
verified that (26) is in agreement with the Kibble-Zurek argument. We also
note that the defect density is overestimated by a factor $\pi$ in [7].

###  2.4 Relation to coarsening process

The correlation function $G_{1}(t)$ displayed in the inset of Fig. 3 for
logarithmic cooling decays fast with time for $Z\gg 1$. But close to the
critical point, it changes very slowly (“frozen”). In view of this behavior,
the Kibble-Zurek argument assumes that the dynamics remain frozen for times
larger than ${t_{*}}$, and the density of defects at the end of quench is
simply inherited from that at $t={t_{*}}=\tau x_{*}$ [2, 20, 23]. However,
from (20) and (21a), we find that

$\frac{D(\tau)}{{\cal D}({t_{*}})}=\frac{1-G_{1}(\tau)}{1-{\cal
G}_{1}({t_{*}})}<1~{},~{}\forall~{}\alpha>0~{}.$ (27)

For logarithmic cooling, the above ratio is $0.59,0.49,0.38$ for
$\alpha=1/2,1,10$, respectively, which shows that the Kibble-Zurek argument
overestimates the density of domain walls by a substantial amount (except, of
course, for $\alpha=0$ for which the remaining time is of order one). It has
been proposed that during the time interval ${t_{*}}<t<\tau$, the dynamics are
not frozen and critical coarsening takes place that results in the decrease in
the number of domain walls [6].

To investigate this proposal, we first note that if the system is
instantaneously quenched from an initial temperature $T_{i}=T({t_{*}})$ to a
final temperature $T_{f}=T(\tau)$, the correlation function is given by [36],

$\displaystyle{\hat{\cal G}}_{k}(t)$ $\displaystyle=$ $\displaystyle
y_{f}^{k}+\frac{2}{\pi}\int_{0}^{\pi}dqe^{-2(1-\gamma\cos(q))(t-{t_{*}})}\sin(kq)$
(28) $\displaystyle\times$
$\displaystyle\left[\frac{y_{i}\sin(q)}{1-2y_{i}\cos(q)+y_{i}^{2}}-\frac{y_{f}\sin(q)}{1-2y_{f}\cos(q)+y_{f}^{2}}\right]~{},~{}t>{t_{*}}~{},$

where $y_{i,f}=\tanh(T_{i,f}^{-1})$. As the above integral gets contribution
from small $q$ and $y_{f}=1$, the first term in the bracket on the RHS of the
above equation can be neglected in comparison to the second term. On carrying
out the resulting integral, we find that

$\displaystyle{\hat{\cal G}}_{k}(t)\approx
1-\frac{k}{\sqrt{\pi(t-{t_{*}})}}~{},~{}~{}k\ll\sqrt{t-{t_{*}}}~{}.$ (29)

For $K=\frac{k}{\sqrt{\tau(1-x_{*})}},Z=\frac{1-x}{1-x_{*}}$ where $1-x_{*}$
depends on the details of the cooling scheme, the above can be written as

${\hat{\cal
G}}_{k}(t)=1-\frac{K}{\sqrt{\pi}}-\frac{KZ}{2\sqrt{\pi}}~{},~{}K,Z\ll 1~{}.$
(30)

Assuming that critical coarsening drives the dynamics for $t>{t_{*}}$ [6],
$G_{k}(t)$ may be approximated by ${\hat{\cal G}}_{k}(t)$. The above equation
then shows that the density of defects at the end of quench has the same
$\tau$-scaling as obtained by assuming that the dynamics remain frozen; this,
however, is a consequence of the fact that the exponent $z={\hat{z}}=2$ in
this model. As expected, the ratio of the density of defects,
$(1-G_{1}(\tau))/(1-{\hat{\cal G}}_{1}(\tau))\to 1$ as $\alpha$ or
$\beta\to\infty$ but, in general, differs from one. Importantly, (30) shows
that close to the critical point, the correlation function ${\hat{\cal
G}}_{k}(t)$ depends linearly on $Z$. Below we therefore study the scaling
behavior of the approach to the critical point to test how well the
instantaneous quench model describes the dynamics when ${t_{*}}<t<\tau$.

####  2.4.1 Logarithmic cooling

For small $K$ and nonzero $Z$, from (17), we obtain

$\displaystyle{\cal G}_{k}(t)$ $\displaystyle-$ $\displaystyle G_{k}(t)=\alpha
K\int_{0}^{\infty}dY(Z+Y)^{\alpha-1}e^{-\frac{2}{\alpha+1}\left[(Z+Y)^{\alpha+1}-Z^{\alpha+1}\right]}$
(31) $\displaystyle\times$
$\displaystyle\left[\frac{1+4Y(Z+Y)^{\alpha}}{\sqrt{2(Z+Y)^{\alpha}}}e^{2Y(Z+Y)^{\alpha}}\textrm{erfc}(\sqrt{2Y(Z+Y)^{\alpha}})-2\sqrt{\frac{Y}{\pi}}\right]~{}.$

On changing the dummy variable from $Y$ to
$w=\frac{2}{\alpha+1}(Z+Y)^{\alpha+1}$ in the above integral and denoting the
resulting integrand by $e^{-w}H(Z,w)$, we obtain

$\displaystyle{\cal G}_{k}(t)-G_{k}(t)$ $\displaystyle=$ $\displaystyle
K\alpha\int_{\frac{2Z^{\alpha+1}}{\alpha+1}}^{\infty}dwe^{-w}H(Z,w)$ (32)
$\displaystyle\stackrel{{\scriptstyle Z\to 0}}{{\approx}}$ $\displaystyle
K\alpha\int_{\frac{2Z^{\alpha+1}}{\alpha+1}}^{\infty}dwe^{-w}\left(H(0,w)+Z\frac{\partial
H(0,w)}{\partial w}\right)~{},$ (33)

where the subleading terms in the above integrand are of order $Z^{1+\alpha}$
and $Z^{2}$ for $\alpha\leq 1$ and $\alpha>1$, respectively. As the integral
over the first term in the parentheses is equal to
$1-G_{k}(\tau)-\sqrt{2}KZ^{\alpha/2}$, we find that to leading order in $Z$,

$\displaystyle G_{k}(\tau)-G_{k}(t)=$ $\displaystyle\frac{K\alpha
Z}{\pi}\int_{\frac{2Z^{\alpha+1}}{\alpha+1}}^{\infty}dw\left(\frac{2}{(\alpha+1)w}\right)^{\frac{3}{2(\alpha+1)}}{I}(w)~{}$
, $\alpha<1/2$ (34a) $\displaystyle G_{k}(\tau)-G_{k}(t)=$
$\displaystyle\frac{K\alpha
Z}{\pi}\int_{0}^{\infty}dw\left(\frac{2}{(\alpha+1)w}\right)^{\frac{3}{2(\alpha+1)}}{I}(w)~{}$
, $\alpha>1/2$ (34b)

where

$I(w)=\sqrt{\pi}e^{-w}(\alpha w+w+1)-\frac{\pi}{2}e^{\alpha
w}(2(\alpha+1)w+3)\sqrt{(\alpha+1)w}\text{erfc}\left(\sqrt{(\alpha+1)w}\right)~{}.$
(35)

On performing the above integrals, we finally arrive at

$\displaystyle G_{k}(\tau)-G_{k}(t)\stackrel{{\scriptstyle Z\ll 1}}{{=}}$
$\displaystyle\frac{4K\alpha}{\sqrt{\pi}(1-2\alpha)}Z^{\alpha+\frac{1}{2}}~{}$
, $\alpha<1/2$ (36a) $\displaystyle
G_{k}(\tau)-G_{k}(t)\stackrel{{\scriptstyle Z\ll 1}}{{=}}$
$\displaystyle\frac{2K\alpha}{\sqrt{\pi}}\ln(1/Z)Z~{}$ , $\alpha=1/2$ (36b)
$\displaystyle G_{k}(\tau)-G_{k}(t)\stackrel{{\scriptstyle Z\ll 1}}{{=}}$
$\displaystyle\frac{K}{2\sqrt{\pi}}\left(\frac{2}{1+\alpha}\right)^{\frac{3}{2(1+\alpha)}}\Gamma\left(\frac{2\alpha-1}{2\alpha+2}\right)Z~{}$
, $\alpha>1/2$ . (36c)

We first verify that (36c) matches (30) when $\alpha\to\infty$ and the
numerical data for finite $\alpha>1/2$ as shown in the inset of Fig. 4. For
$\alpha<1/2$, to see if (36a) agrees with numerical results, we need to
consider very large $\tau$ and very small $Z$ since the subleading corrections
to $G_{1}(\tau)$ (as already discussed) and to the integral in (34a) (see
below (33)) decay slowly for small $\alpha$. Figure 4 shows that in accordance
with (36a), $(G_{1}(\tau)-G_{1}(t))/Z$ diverges as $Z\to 0$ and approaches the
prefactor with increasing $\tau$. On comparing (30) with (36a)-(36c), we
arrive at the conclusion that the instantaneous quench model captures the
dependence on $Z$ quantitatively for $\alpha\geq 1/2$.

####  2.4.2 Algebraic cooling

For large $\tau$, the equilibrium correlation function ${\cal
G}_{k}(t)=e^{-k\sqrt{2(1-\gamma)}}$ is equal to one for small $Z$. Comparing
(25a) and (30), we find that for any $\beta$, the defect density is exactly
given by the instantaneous quench model; this is because of the fact that
algebraic cooling corresponds to $\alpha\to\infty$ in logarithmic cooling.

##  3 Zero-range process in mean-field geometry

In this section, we consider a zero-range process in mean-field geometry where
each site is connected to every other site. A site can have $m\geq 0$
identical particles each with mass one, and a particle can hop out to another
site with a rate $u(m,t)$ that depends on the number of particles present at
the departure site at time $t$. The single-site mass distribution $P(m,t)$ at
time $t$ evolves according to the following equation,

$\displaystyle\frac{dP(0,t)}{dt}$ $\displaystyle=$ $\displaystyle
u(1,t)P(1,t)-w(t)P(0,t)~{},$ (37) $\displaystyle\frac{dP(m,t)}{dt}$
$\displaystyle=$ $\displaystyle u(m+1,t)P(m+1,t)+w(t)P(m-1,t)$ (38)
$\displaystyle-$ $\displaystyle[u(m,t)+w(t)]P(m,t)~{},~{}m\geq 1~{},$

where $w(t)=\sum_{m^{\prime}=1}^{\infty}u(m^{\prime},t)P(m^{\prime},t)$ is the
mean hop rate. In the following, we will work with the following time-
dependent hop rate,

$u(m,t)=1+\frac{b(t)}{m}~{},~{}m>0~{},$ (39)

where

$b(t)=b_{c}\left[1-\left(1-x\right)^{\alpha}\right]~{},~{}x=\frac{t}{\tau}$
(40)

that changes from zero to a final value $b_{c}$ algebraically in time at rate
$\tau^{-1}$. In the above equation, $b_{c}>2$ (see below) and, as before, the
exponent $\alpha>0$. In this section, we will not consider the annealing
scheme analogous to (9b) since it corresponds to the instantaneous quench
limit $\alpha\to\infty$, as seen in the last section. We also note that unlike
in the Ising model studied in the last section where the correlation function
obeys (5), here the equations of motion for the mass distribution are
nonlinear (due to the fugacity term) and the coefficients are also mass-
dependent.

When the hop rates are time-independent and given by $u(m)=1+(b/m)$, a
stationary state exists and exhibits a phase transition from a fluid phase in
which particles are homogeneously distributed to a condensate phase where a
finite fraction of particles reside in a single mass cluster, as the parameter
$b$ is increased keeping the total particle density $\rho_{c}$ constant [34,
9]. The stationary state mass distribution is known exactly to be

${\cal P}(m)=\frac{\omega^{m}f(m)}{g(\omega)}~{},$ (41)

where $f(m)=\prod_{k=1}^{m}u^{-1}(k)$ and the partition function,
$g(\omega)=\sum_{m=0}^{\infty}\omega^{m}f(m)={{}_{2}}F_{1}(1,1;1+b;\omega)$ is
the Gauss hypergeometric function. The total mass density $\rho_{c}$ is
related to the fugacity $\omega=\sum_{m=1}^{\infty}u(m){\cal P}(m)$ as
$\rho_{c}=\omega g^{\prime}(\omega)/g(\omega)$ where prime denotes derivative
with respect to $\omega$. For the above $u(m)$, with increasing $b$, the
fugacity reaches its maximum value one at a critical value,
$b_{c}=2+\rho_{c}^{-1}~{},~{}b_{c}>2$. Close to the critical point, the steady
state distribution ${\cal P}(m)\sim m^{-b}e^{-m/\xi}$ where the static
correlation length, $\xi\sim(1-\omega)^{-1}\sim(b_{c}-b)^{-\nu}$ with [39, 9]

$\displaystyle{\nu}=$ $\displaystyle(b_{c}-2)^{-1}~{}$ , $2<b_{c}<3$ (42a)
$\displaystyle{\nu}=$ $\displaystyle 1~{}$ , $b_{c}\geq 3$ . (42b)

For later reference, we also note that for this model, the stationary state
dynamical exponent $z=2$ at the critical point [40] and the critical
coarsening exponent ${\hat{z}}=2$ [35].

In the following subsections, we quench the system from $b(0)=0$ to
$b(\tau)=b_{c}$ at a fixed density $\rho_{c}$, and study how the mass
distribution behaves when the hop rates are time-dependent. For convenience,
we will work with the ratio $R(m,t)\equiv P(m,t)/{\cal P}(m,t)$ of
nonequilibrium to equilibrium probability distribution (instead of the
difference between them).

###  3.1 Dynamics deep in the fluid phase

Our numerical solution of (37) and (38) suggests that the mass distribution
$P(m,t)$ has different scaling form for small and large mass. Therefore, for
large $t,\tau$ with $x=t/\tau$ finite, we make the following scaling ansatz,

$\displaystyle\frac{P(m,t)}{{\cal P}(m,t)}=$ $\displaystyle
1+\tau^{-\theta}\delta P(m,x)~{}$ ,
$m\epsilon\stackrel{{\scriptstyle\tau\to\infty}}{{\to 0}}$ (43a)
$\displaystyle\frac{P(m,t)}{{\cal P}(m,t)}=$ $\displaystyle Q(M,x)~{}$ ,
$m\epsilon$ finite , (43b)

where the scaling variable $M=m\epsilon(\tau)$ for large $m$ and
$\epsilon(\tau)$ decays as a power law with $\tau$. For consistency with
(43a), we also have $Q(0,x)=1$. Using the above scaling ansatz in the
definition of the fugacity, we find that

$\displaystyle w(t,\tau)$ $\displaystyle=$
$\displaystyle\int_{0}^{\epsilon^{-1}}dm{\cal P}(m,x)(1+\tau^{-\theta}\delta
P(m,x))u(m,x)$ (44) $\displaystyle+$
$\displaystyle\epsilon^{b_{c}-1}\int_{1}^{\infty}dMM^{-b_{c}}e^{-\frac{M}{\epsilon\xi(x)}}Q(M,x)$
$\displaystyle\approx$ $\displaystyle{\omega(x)}[1+\tau^{-\theta}\delta
w(x)]~{},$ (45)

since the second integral on the RHS of (44) vanishes in the $\tau\to\infty$
limit as this term is exponentially small in $\epsilon^{-1}$. Thus the
contribution to the fugacity comes from small mass, as one may expect in the
fluid phase which is characterized by typical mass of order unity.

Using (43a) and (45) in the master equation (38) for $m\ll\epsilon^{-1}$, we
find that the terms on the left-hand side (LHS) are of order $\tau^{-1}$ while
the RHS is of order $\tau^{-\theta}$. We therefore deduce that the exponent
$\theta=1$ which is in line with the data collapse shown in the inset of Fig.
5 for fugacity. For large mass, using the scaling ansatz (43b) in the master
equation (38) and taking the scaling limits $m\to\infty,\epsilon\to 0$ with
finite $M$, we find that the scaling function $Q(M,x)$ obeys the following
differential equation,

$\displaystyle\frac{MQ(M,x)}{\tau\epsilon(\tau)}\frac{d\ln\omega(x)}{dx}=\epsilon(\tau)(\omega(x)-1)\frac{\partial
Q(M,x)}{\partial M}~{},$ (46)

where we have also used the exact stationary state distribution (41). The
above equation yields $\epsilon\sim\tau^{-1/2}$ and the scaling function

$Q(M,x)=\exp\left[-\frac{m^{2}}{2\tau\omega(x)(1-\omega(x))}\frac{d\omega(x)}{dx}\right]~{}.$
(47)

Figure 5 shows that the numerical solution of (37) and (38) for the mass
distribution in the fluid phase is in good agreement with (47).

###  3.2 Dynamics close to the critical point

For $x\to 1$, we make the following scaling ansatz,

$\displaystyle\frac{P(m,t)}{{\cal P}(m,t)}=$ $\displaystyle
1+\tau^{-\theta}\delta P(m,Z)~{}$ ,
$m\epsilon\stackrel{{\scriptstyle\tau\to\infty}}{{\to}}0$ (48a)
$\displaystyle\frac{P(m,t)}{{\cal P}(m,t)}=$ $\displaystyle Q(M,Z)~{}$ ,
$m\epsilon$ finite , (48b)

where, besides the scaled mass $M=m\epsilon(\tau)$, we have the scaled
remaining time $Z=(1-x)\Lambda(\tau)$, $\Lambda(\tau)$ being an algebraically
increasing function of $\tau$. Although the correlation length diverges close
to the critical point, as we shall see below, $\epsilon\xi$ remains finite and
$\tau^{-\theta}\leq\epsilon^{b_{c}-1}$. Then the second term on the RHS of
(44) for fugacity vanishes as $\epsilon\to 0$ since $b_{c}>2$ (assuming the
integral is finite). It thus follows that for large $\tau$, the fugacity
$w(Z)=\omega(x)[1+\tau^{-\theta}\delta w(Z)]$, where the exponent $\theta$ is
determined below.

For large mass, we use the scaling ansatz (48b) in the master equation (38)
and find that

$\displaystyle\frac{\epsilon^{2}}{2}\left(\omega+\frac{wu(m,x)}{\omega}\right)\frac{\partial^{2}Q}{\partial
M^{2}}+\epsilon\left(\omega-\frac{wu(m,x)}{\omega}\right)\frac{\partial
Q}{\partial M}$ (49) $\displaystyle+$
$\displaystyle\left(\omega-u(m,x)-w+\frac{wu(m,x)}{\omega}\right)Q=\frac{1}{\tau}\frac{d\ln{\cal
P}(m,x)}{dx}Q-\frac{\Lambda}{\tau}\frac{\partial Q}{\partial Z}~{}.$

Close to the critical point, we also have

$\displaystyle\omega-u(m,x)$ $\displaystyle\approx$
$\displaystyle\omega-1-\frac{b_{c}}{m}$ (50) $\displaystyle 1-\omega(x)$
$\displaystyle\approx$ $\displaystyle C(b_{c})(1-x)^{\alpha\nu}$ (51)
$\displaystyle\frac{d\ln{\cal P}(m,x)}{dx}$ $\displaystyle\approx$
$\displaystyle(m-\rho_{c})C\alpha\nu(1-x)^{\alpha\nu-1}+\alpha
b_{c}(H_{b_{c}}-H_{b_{c}+m})(1-x)^{\alpha-1}~{},$ (52)

where $H_{n}$ is harmonic number. Using these stationary state properties in
(49), a simple power counting yields
$\Lambda\sim\tau^{\frac{1}{1+2\alpha\nu}}$ and
$\epsilon\sim\Lambda^{-\alpha\nu}$ so that we finally have

$\displaystyle\frac{P(m,t,\tau)}{{\cal P}(m,x)}=$ $\displaystyle
Q(m\tau^{-\frac{\alpha}{b_{c}+2\alpha-2}},(1-x)\tau^{\frac{b_{c}-2}{b_{c}+2\alpha-2}})~{}$
, $2<b_{c}<3$ (53a) $\displaystyle\frac{P(m,t,\tau)}{{\cal P}(m,x)}=$
$\displaystyle
Q(m\tau^{-\frac{\alpha}{1+2\alpha}},(1-x)\tau^{\frac{1}{1+2\alpha}})~{}$ ,
$b_{c}\geq 3$ . (53b)

The above scaling laws are consistent with the Kibble-Zurek argument which
predicts $\xi\sim\tau^{\frac{\alpha\nu}{1+\alpha\nu
z}},1-x\sim\tau^{-\frac{1}{1+\alpha\nu z}}$ since $z=2$ (see also Fig. 6). (We
also verify the above claim that $\epsilon\xi\sim\epsilon(1-\omega)^{-1}$ is
of order one.)

The above analysis also shows that the scaling function $Q(M,Z)$ obeys the
following differential equation,

$\frac{\partial^{2}Q}{\partial
M^{2}}-\left[CZ^{\alpha\nu}+\frac{b_{c}}{M}+\frac{\tau^{-\theta}\delta
w}{\epsilon}\right]\frac{\partial Q}{\partial M}+\frac{\tau^{-\theta}\delta
w}{\epsilon}\left[CZ^{\alpha\nu}+\frac{b_{c}}{M}\right]Q=-\frac{\partial
Q}{\partial Z}+\alpha\nu CMZ^{\alpha\nu-1}Q~{},$ (54)

where the constant $C$ is given by (see Appendix A),

$\displaystyle C(b_{c})=$ $\displaystyle\left(\frac{b_{c}\sin(\pi
b_{c})}{\pi(b_{c}-2)^{2}(b_{c}-1)^{2}}\right)^{\frac{1}{b_{c}-2}}~{}$ ,
$2<b_{c}<3$ (55a) $\displaystyle C(b_{c})=$
$\displaystyle\frac{b_{c}(b_{c}-3)}{(b_{c}-1)^{2}}~{}$ , $b_{c}>3$ . (55b)

The $\delta w$ term in the above differential equation contributes if
$\tau^{-\theta}\sim\epsilon$ which, as shown below in (57b), is true for
$b_{c}<3$. Thus we obtain a closed equation for $Q(M,Z)$ for $b_{c}>3$ whose
numerical solution subject to boundary conditions
$Q(0,Z)=1,Q(\infty,Z)=0,Q(M,\infty)=0$ matches well with that obtained using
the master equation, as shown in Fig. 6.

To find the exponent $\theta$, we adapt the analysis of [35] for the
coarsening dynamics of zero-range process with time-independent hop rates to
the model considered here. Using (48a) in the master equation (38) for small
mass, we arrive at

$\Delta P(m+1,Z)-u(m,1)\Delta P(m,Z)+\delta
w(x)(u(m,1)-1)\approx\frac{\tau^{\theta}}{\tau}\frac{d\ln{\cal
P}(m,x)}{dx}~{},$ (56)

where $\Delta P(m,Z)=\delta P(m,Z)-\delta P(m-1,Z)$. While the LHS is of order
unity, the RHS scales as $m\tau^{\theta-1}(1-x)^{\alpha-1}\sim
m\tau^{\theta-\frac{\alpha(1+2\nu)}{1+2\alpha\nu}}$ (see (52)). Applying (56)
to small enough mass so that the RHS vanishes for large $\tau$, we find that
$\delta P(m,Z)=m\delta w(Z)+\delta P(0,Z)$. Furthermore, using (48a) and (48b)
in the total mass density equation, $\rho_{c}=\int_{0}^{\infty}dmP(m,t)$, we
find that $\tau^{-\theta}\int_{0}^{\epsilon^{-1}}dmm^{2}{\cal
P}(m,x)\sim\epsilon^{b_{c}-2}$. As the stationary state mass variance close to
the critical point is finite for $b_{c}>3$ and diverges as
$\epsilon^{b_{c}-3}$ for $2<b_{c}<3$, we find that the fugacity exponent,

$\displaystyle\theta=$ $\displaystyle\frac{\alpha}{b_{c}-2+2\alpha}~{}$ ,
$2<b_{c}<3$ (57a) $\displaystyle\theta=$
$\displaystyle\frac{\alpha(b_{c}-2)}{1+2\alpha}~{}$ , $b_{c}\geq 3$ , (57b)

on using (42a) and (42b). We verify that the result (27) of [35] for
instantaneous quench to the critical point when $b_{c}>3$ is recovered when
$\alpha\to\infty$; moreover, the above equation predicts the corresponding
exponent for $2<b_{c}<3$ to be $1/2$. The exponent $\theta$ is numerically
tested in Fig. 7, and we find a good agreement with (57a) and (57b).

###  3.3 Relation to coarsening process

To understand the relationship to coarsening process, using (28) of [35], we
first note that

$\displaystyle\frac{{\hat{\cal P}}(m,t)}{{\cal P}(m,t)}$ $\displaystyle=$
$\displaystyle
1-\frac{1}{2^{1+b_{c}}\Gamma(\frac{3+b_{c}}{2})}\left(\frac{m}{\sqrt{t-t_{*}}}\right)^{1+b_{c}}$
(58) $\displaystyle=$ $\displaystyle
1-\frac{1}{2^{1+b_{c}}\Gamma(\frac{3+b_{c}}{2})}\left(\frac{M}{\sqrt{1-Z}}\right)^{1+b_{c}}$
(59) $\displaystyle\approx$ $\displaystyle
1-\frac{M^{1+b_{c}}}{2^{1+b_{c}}\Gamma(\frac{3+b_{c}}{2})}-\frac{M^{1+b_{c}}Z}{2^{1+b_{c}}\Gamma(\frac{1+b_{c}}{2})}~{},$
(60)

where ${\hat{\cal P}}(m,t)$ denotes the mass distribution on instantaneous
quench. In other words, close to the critical point, the ratio of the mass
distribution decays linearly with the scaled remaining time $Z$, when the
system is quenched infinitely fast.

To see if this linear behavior holds when quenching occurs at a finite rate,
we studied the ratio $R(m,t)=P(m,t)/{\cal P}(m,t)$ numerically. As our results
displayed in Fig. 8 show, $R(m,t)$ decays linearly with $Z$ for $\alpha=1$ but
not for $\alpha=1/2$. Although we are not able to show it analytically, our
numerical results suggest that the linear scaling holds for $\alpha\geq 1$ and
$R(m,t)-R(m,\tau)$ varies as $Z^{\alpha}$ for $\alpha\leq 1$. We have also
studied how the mass variance changes with $Z$ and find that it also varies
sublinearly for sufficiently small $\alpha$ (data not shown).

##  4 Discussion

In this article, we have analyzed in detail the nonequilibrium dynamics of a
classical system when it is annealed slowly from a disordered phase to the
critical point in the framework of a kinetic Ising model and a zero-range
process. The Kibble-Zurek argument that explains how the equilibrium is
approached with decreasing annealing rate has been verified numerically in
various recent studies [8, 26, 9, 10, 11]. But it has also been found that
this argument overestimates the defect density, and scaling laws different
from those predicted by it can be obtained when critical coarsening is taken
into account [6]. However, the previous body of work does not give any insight
into the associated scaling functions.

Here, using the dynamical equations for the correlation function in the Ising
model and mass distribution in the zero-range process, we have derived these
scaling laws and find that Kibble-Zurek scaling holds; we also find that
scaling function obeys (B.3) for the Ising model and (54) for the zero-range
process. To investigate the role of critical coarsening in the slow quench
dynamics, we have focused on the scaling of appropriate quantities in these
models with the remaining time $Z$. Working with the annealing scheme (9a) and
(40), we find that close to the critical point, ${\cal V}$ which represents
domain wall density in the Ising model and nonequilibrium mass distribution
(relative to the equilibrium one) in zero-range process is of the following
form,

$\displaystyle{\cal V}(Z)=$ $\displaystyle
A_{1}(\tau,\alpha)+A_{2}(\tau,\alpha)Z^{a}~{},$ $\alpha<\alpha_{c}$ (61a)
$\displaystyle{\cal V}(Z)=$
$\displaystyle{\hat{A}}_{1}(\tau,\alpha)+{\hat{A}}_{2}(\tau,\alpha)Z~{},$
$\alpha>\alpha_{c}$ , (61b)

where the exponent $a<1$, and the critical quench exponent $\alpha_{c}=1/2$
and $1$ in the Ising model and zero-range process, as depicted in Fig. 4 and
Fig. 8, respectively. As the linear scaling is obtained following rapid quench
(see (30) and (60)), we conclude that critical coarsening plays a role when
the final quench point is approached sufficiently fast.

We close this article with some open questions and directions. It would be
instructive to study the slow quench dynamics for the Ising chain with
Kawasaki dynamics where, unlike for Glauber dynamics studied here, the dynamic
exponent $z$ and coarsening exponent ${\hat{z}}$ are not equal, and one may
expect the defect density at the end of quench to exhibit a scaling different
from that predicted by Kibble-Zurek argument [6]. Here we have quenched the
system to the critical point, but a nontrivial critical point $b_{c}$ exists
in the zero-range process and a detailed understanding of slow quench dynamics
in the ordered phase will be carried out in a future study. More general
schemes such as series of cooling-heating cycles that in the context of
glasses [4] are known to show hysteretic effect could also be interesting to
study. Most of the previous work on high-dimensional models have assumed
linear annealing [6, 9]; it would be interesting to explore the dependence on
other functional forms and determine the critical quench exponent above which
instantaneous quench model works. An analytical understanding of the mean-
field zero-range process presented here, and more generally, high-dimensional
models with nontrivial critical point remains an open question for future
work.

##  Appendix A Zero-range process: density-fugacity relation

From the density-fugacity relation $\rho=\omega g^{\prime}(\omega)/g(\omega)$,
we have

$\rho_{c}=\frac{1}{b_{c}-2}=\frac{\omega}{b+1}\frac{{{}_{2}}F_{1}(2,2;2+b;\omega)}{{{}_{2}}F_{1}(1,1;1+b;\omega)}~{}.$
(A.1)

Using (15.3.6) of [41], we expand the RHS of the above equation about
$\omega=1$. For $b>3$, we get

$\frac{1}{b_{c}-2}\approx\frac{1}{b-2}\left(1-\frac{(b-1)^{2}(1-\omega)}{(b-2)(b-3)}\right)$
(A.2)

which yields (55b) on using the annealing scheme (40) when $x\to 1$.
Similarly, for $2<b<3$, we have

$\frac{1}{b_{c}-2}\approx\frac{1}{b-2}-\frac{\pi(b-1)^{2}}{\sin(\pi
b)}(1-\omega)^{b-2}$ (A.3)

which leads to (55a).

##  Appendix B Ising model: scaling functions

For a comparison with the differential equations (46) and (54) for the scaling
functions, respectively, in the fluid phase and close to the critical point in
the zero-range process, here we give the corresponding differential equations
for the Ising model. In the paramagnetic phase, the ratio of the correlation
functions $R(K,x)=\frac{G_{k}(t)}{{\cal G}_{k}(t)}$ obeys

$\frac{R(K,x)}{\tau}\frac{\partial\ln{\cal C}(k,x)}{\partial
x}=-2\sqrt{1-\gamma^{2}}\epsilon(\tau)\frac{\partial R(K,x)}{\partial K}~{},$
(B.1)

which leads to

$\displaystyle\frac{G_{k}(t)}{{\cal
G}_{k}(t)}=\exp\left[-\frac{k^{2}}{\tau}\frac{\frac{d\gamma(x)}{dx}}{4\gamma(x)(1-\gamma^{2}(x))}\right]~{}.$
(B.2)

Close to the critical point, the ratio $R(K,Z)=\frac{G_{k}(t)}{{\cal
G}_{k}(t)}$ obeys the following equation,

$\frac{\partial^{2}R}{\partial K^{2}}-2\sqrt{2}Z^{\alpha/2}\frac{\partial
R}{\partial K}=-\frac{\partial R}{\partial
Z}+\frac{K\alpha}{\sqrt{2}}Z^{\frac{\alpha}{2}-1}R~{},$ (B.3)

with $R(0,Z)=1,R(\infty,Z)=0,R(K,\infty)=1$ for logarithmic cooling.
Alternatively, the scaling function $F(K,Z)=(1-R(K,Z)){\cal G}(K,Z)$ obeys the
following equation,

$2Z^{\alpha}F-\frac{\partial F}{\partial Z}-\frac{\partial^{2}F}{\partial
K^{2}}=\frac{\alpha}{\sqrt{2}}KZ^{\frac{\alpha}{2}-1}{\cal G}(K,Z)~{}.$ (B.4)

We have checked that (17) satisfies (B.4).

Figure 1: Kinetic Ising model: Deviation in the correlation function, ${\cal
G}_{k}(t)-G_{k}(t)$ for logarithmic cooling scheme (9a) with $\tau=2\times
10^{3}$ and $\alpha=3/2$ at
$t=500(\vartriangle),1000(\square),1500(\Circle),1800(\triangledown),1900(\lozenge),2000(\circ)$
obtained by numerically solving the dynamical equation (5) and using the
equilibrium correlation function (6). The lines show the expression (13) for
$t\leq 1500$ in the main figure. The deviation ${\cal G}_{k}(t)-G_{k}(t)$
increases with time, and is substantial when the system is close to the
critical point (see inset).

Figure 2: Kinetic Ising model: Nonequilibrium correlation function,
$G_{k}(\tau)$ at the critical point for logarithmic cooling scheme (9a) and
$\alpha=3/2$ with $\tau=2\times 10^{3}(\vartriangle),10^{4}(\square),5\times
10^{4}(\circ)$. The points show the data obtained by numerically solving the
dynamical equation (5), and the line shows the expression (21a) for small $k$.

Figure 3: Kinetic Ising model: Defect density $D(\tau)$ at the critical point
for logarithmic cooling scheme (9a) and $\tau=2000$ as a function of the
quench exponent $\alpha$. The numerical data ($\vartriangle$) does not agree
well with the analytical result (21a) ($\triangledown$) for small $\alpha$
since the correction $\delta D(\tau)$ to the leading order result decays
slowly (see bottom, left inset for $\alpha=1/4$). The top, right inset shows
the correlation function $G_{1}(t)$ for $\tau=2000$ and $\alpha=1/4,1/2,3/4$
(bottom to top) as a function of the scaled time
$Z=(1-x)\tau^{\frac{1}{1+\alpha}}$. Note that $G_{1}(t)$ changes very slowly
for $Z<1$.

Figure 4: Kinetic Ising model: (Scaled) correlation function
${G}_{1}(\tau)-G_{1}(t)$ close to the critical point as a function of
$Z=(1-x)\tau^{\frac{1}{1+\alpha}}$ for logarithmic cooling scheme (9a) and
$\alpha=1/4$ with $\tau=10^{8}$ (dotted) and $10^{10}$ (solid) obtained by
solving (5) numerically. The dashed curve shows the analytical result (36a).
As the curves diverge for $Z\to 0$, it follows that the correlation function
$G_{1}(t)$ decays sublinearly with $Z$ for $\alpha<1/2$. In contrast, the
nearest neighbor correlation exhibits linear decay for $\alpha>1/2$, as the
inset shows for $\alpha=3/2$ with $\tau=10^{4}$ for numerical data (dotted)
and analytical expression (36c) (dashed).

Figure 5: Zero-range process: Ratio of the mass distribution function,
$R(m,t)=P(m,t)/{\cal P}(m,t)$ for large mass in the fluid phase as a function
of scaled mass $M=m/\sqrt{\tau}$ for annealing scheme (40) with $\alpha=1$.
The points show the data obtained by numerically solving the master equation
(37) and (38) for $x=1/2$, and the solid line shows the scaling function
$Q(M,x)$ given by (47). The inset figure shows the data collapse for scaled
fugacity difference $\delta w$ far from the critical point to support the
conclusion that the exponent $\theta=1$ in (45). In both the figures,
$\tau=2^{16}(\circ),2^{17}(\vartriangle),2^{18}(\Box)$ and $b_{c}=5/2$.

Figure 6: Zero-range process: Data collapse for the ratio of the mass
distribution functions, $R(m,t)=P(m,t)/{\cal P}(m,t)$ for large mass and close
to the critical point to confirm the scaling (53a) and (53b) for annealing
scheme (40) with $\alpha=1$. The remaining time $Z\approx 0.112$ for $b_{c}=4$
(main) and $Z\approx 0.051$ for $b_{c}=5/2$ (inset), and
$\tau=2^{17}(\circ),2^{18}(\vartriangle),2^{19}(\Box)$. The solid line in the
main figure shows the numerical solution of the differential equation (54) for
the scaling function.

Figure 7: Zero-range process: Data collapse for the fugacity difference,
$\delta w(Z)=\tau^{\theta}(w-\omega)$ for annealing scheme (40) with
$\alpha=1$, $b_{c}=5/2$ (main) and $7/2$ (inset), and
$\tau=2^{17}(\circ),2^{18}(\vartriangle),2^{19}(\Box)$ to verify the
$\theta$-exponent given by (57a) and (57b).

Figure 8: Zero-range process: Variation of the ratio of the mass distribution,
$R(m,t)=P(m,t)/{\cal P}(m,t)$ close to the critical point with the scaled
remaining time $Z$ for $M=2,b_{c}=5/2(\vartriangle)$ and
$M=4,b_{c}=7/2(\circ)$ with $\tau=2^{20}$ for annealing scheme (40). The
difference $R(m,t)-R(m,\tau)$ decays sublinearly with $Z$ for $\alpha=1/2$
(main), and linearly for $\alpha=1$ (inset) as indicated by solid line and in
agreement with the instantaneous quench model prediction (60).

## References

  * [1] Bray A J 1994 Adv. Phys. 43 357
  * [2] Kibble T W B 1976 J. Phys. A 9 1387
  * [3] Cornell S J, Kaski K and Stinchcombe R B 1991 Phys. Rev. B 44 12263
  * [4] Ritort F and Sollich P 2003 Adv. Phys. 52 219–342
  * [5] Cornell S J, Kaski K and Stinchcombe R 1992 Phys. Rev. B 45 2725–2738
  * [6] Biroli G, Cugliandolo L F and Sicilia A 2010 Phys. Rev. E 81 050101(R)
  * [7] Krapivsky P L 2010 J. Stat. Mech. P02014
  * [8] Jelić A and Cugliandolo L 2011 J. Stat. Mech. P02032
  * [9] Priyanka and Jain K 2016 EPL 116 26003
  * [10] Ricateau H, Cugliandolo L F and Picco M 2018 J. Stat. Mech. P013201
  * [11] Jeong K, Kim B and Lee S J 2019 Phys. Rev. E 99 022113
  * [12] Mathey S and Diehl S 2020 Phys. Rev. Research 2 013150
  * [13] Su Z and Millis A J 2020 Phys. Rev. X 10 021028
  * [14] Zurek W H, Dorner U and Zoller P 2005 Phys. Rev. Lett. 95 105701
  * [15] Polkovnikov A 2005 Phys. Rev. B 72 161201(R)
  * [16] Mukherjee V, Divakaran U, Dutta A and Sen D 2007 Phy. Rev. B 76 174403
  * [17] Dziarmaga J 2010 Adv. Phys. 59 1063–1189
  * [18] Chandran A, Erez A, Gubser A and Sondhi A 2012 Phys. Rev. B 86 064304
  * [19] Karevski D and Harris R J 2016 J. Stat. Mech.: Theo. and Exp. 2016 033204
  * [20] Zurek W H 1985 Nature 317 505
  * [21] del Campo A and Zurek W H 2014 Int. J. Mod. Phys. A 29 1430018
  * [22] Beugnon J and Navon N 2017 J. Phys. B: At. Mol. Opt. Phys. 50 022002
  * [23] Zurek W H 1996 Phys. Rep. 276 177
  * [24] Laguna P and Zurek W H 1997 Phys. Rev. Lett. 78 2519
  * [25] Yates A and Zurek W H 1998 Phys. Rev. Lett. 80 5477
  * [26] Liu C W, Polkovnikov A and Sandvik A W 2014 Phys. Rev. B 89 054307
  * [27] Huse D A and Fisher D S 1986 Phys. Rev. Lett. 57 2203
  * [28] Jäckle J, Stinchcombe R B and Cornell S 1991 J. Stat. Phys. 62 425–433
  * [29] Suzuki S 2009 J. Stat. Mech. P03032
  * [30] Xu N, Wu K H, Rubin S J, Kao Y J and Sandvik A W 2017 Phys. Rev. E 96 052102
  * [31] Reiss H 1980 Chemical Physics 47 15
  * [32] Schilling R 1988 J. Stat. Phys. 53 1227
  * [33] Brey J J and Prados A 1994 Phys. Rev. B 49 984–997
  * [34] Evans M R and Hanney T 2005 J. Phys. A: Math. Gen. 38 195
  * [35] Godréche C 2003 J. Phys. A: Math. Gen. 36 6313–6328
  * [36] Glauber R J 1963 J. Math. Phys. 4 294–307
  * [37] Aliev M A 2009 J. Math. Phys. 50 083302
  * [38] Privman V 2005 Nonequilibrium Statistical Mechanics in One Dimension (Cambridge University Press)
  * [39] Priyanka, Ayyer A and Jain K 2014 Phys. Rev. E 90 062104
  * [40] Godréche C and Luck J M 2001 Eur. Phys. J. B 23 473–486
  * [41] Abramowitz M and Stegun I 1964 Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables (Dover)


# Finnish Paraphrase Corpus

Jenna Kanerva, Filip Ginter, Li-Hsin Chang, Iiro Rastas, Valtteri Skantsi,
Jemina Kilpeläinen, Hanna-Mari Kupari, Jenna Saarni, Maija Sevón, and Otto
Tarkka
TurkuNLP group
Department of Computing
Faculty of Technology
University of Turku, Finland
<EMAIL_ADDRESS>
###### Abstract

In this paper, we introduce the first fully manually annotated paraphrase
corpus for Finnish containing 53,572 paraphrase pairs harvested from
alternative subtitles and news headings. Out of all paraphrase pairs in our
corpus 98% are manually classified to be paraphrases at least in their given
context, if not in all contexts. Additionally, we establish a manual candidate
selection method and demonstrate its feasibility in high quality paraphrase
selection in terms of both cost and quality.

##  1 Introduction

The powerful language models that have recently become available in NLP have
also resulted in a distinct shift towards more meaning-oriented tasks for
model fine-tuning and evaluation. The most typical example is entailment
detection, with the paraphrase task raising in interest recently. Paraphrases,
texts that express the same meaning with differing words Bhagat and Hovy
(2013), are — already by their very definition — a suitable target to induce
and evaluate models’ ability to represent meaning. Paraphrase detection and
generation has numerous direct applications in NLP Madnani and Dorr (2010),
among others in question answering Soni and Roberts (2019), plagiarism
detection Altheneyan and Menai (2019), and machine translation Mehdizadeh
Seraj et al. (2015).

Research in paraphrase naturally depends on the availability of datasets for
the task. We will review these in more detail in Section 2, nevertheless,
barring few exceptions, paraphrase corpora are typically large and gathered
automatically using one of several possible heuristics. Typically a
comparatively small section of the corpus is manually classified to serve as a
test set for method development. The heuristics used to gather and filter the
corpora naturally introduce a bias to the corpora which, as we will show later
in this paper, demonstrates itself as a tendency towards short examples with a
relatively high lexical overlap. Addressing this bias to the extent possible,
and providing a corpus with longer, lexically more diverse paraphrases is one
of the motivations for our work. The other motivation is to cater for the
needs of Finnish NLP, and improve the availability of high-quality, manually
annotated paraphrase data specifically for the Finnish language.

In this paper, we therefore aim for the following contributions: Firstly, we
establish and test a fully manual procedure for paraphrase candidate selection
with the aim of avoiding a selection bias towards short, lexically overlapping
candidates. Secondly, we release the first fully manually annotated paraphrase
corpus of Finnish, sufficiently large for model training. The number of
manually annotated examples makes the released dataset one of the largest, if
not the largest manually annotated paraphrase corpus for any language. And
thirdly, we report the experiences, tools, and baseline results on this new
dataset, hopefully allowing other language NLP communities to assess the
potential of developing a similar corpus for other languages.

##  2 Related Work

Corpus | Data source | Size autom. | Size manual | Labels
---|---|---|---|---
English |  |  |  |
MRPC | Online news | — | 5,801 | 0/1
TUC | News tweets | — | 52K | 0/1
ParaSCI | Scientific papers | 350K | — | 1-5
PARADE | flashcards (computer sci.) | — | 10K | 0-3
QQP | Quora | 404K | — | 0/1
Finnish |  |  |  |
Opusparcus | OpenSubtitles | 480K* | 3,703 | 1-4
TaPaCo | Tatoeba crowdsourcing | 12K | — | —
Table 1: Summary of available paraphrase corpora of naturally occurring
sentential paraphrases. The corpora sizes include the total amount of pairs in
the corpus (i.e. also those labeled as non-paraphrases), thus the actual
number of good paraphrases depend on the class distribution of each corpus.
*The highest quality cutpoint estimated by the authors.

Statistics of the different paraphrase corpora most relevant to our work are
summarized in Table 1. For English, the Microsoft Research Paraphrase Corpus
(MRPC) Dolan and Brockett (2005) is extracted from an online news collection
by applying heuristics to recognize candidate document pairs and candidate
sentences from the documents. Paraphrase candidates are subsequently filtered
using a classifier, before the final manual binary annotation (paraphrase or
not). In the Twitter URL Corpus (TUC) Lan et al. (2017), paraphrase candidates
are identified by recognizing shared URLs in news related tweets. All
candidates are manually binary-labeled. ParaSCI Dong et al. (2021) is created
by collecting paraphrase candidates from ACL and arXiv papers using heuristics
based on term definitions, citation information as well as sentence embedding
similarity. The extracted candidates are automatically filtered, but no
manually annotated data is available. PARADE He et al. (2020) is created by
collecting online user-generated flashcards for computer science related
concepts. All definitions for the same term are first clustered, and
paraphrase candidates are extracted only among a cluster to reduce noise in
candidate selection. All extracted candidates are manually annotated using a
scheme with four labels. Quora Question Pairs (QQP)111data.quora.com/First-
Quora-Dataset-\Release-Question-Pairs contains question headings from the
forum with binary labels into duplicate-or-not questions. The QQP dataset is
larger than other datasets, however, although including human-produced labels,
the labeling is not originally designed for paraphrasing and the dataset
providers warn about labeling not guaranteed to be perfect.

Another common approach for automatic paraphrase identification is through
language pivoting using multilingual parallel datasets. Here sentence
alignments are used to recognize whether two different surface realizations
share an identical or near-identical translation, assuming that the identical
translation likely implies a paraphrase. There are two different multilingual
paraphrase datasets automatically extracted using language pivoting,
Opusparcus Creutz (2018) and TaPaCo Scherrer (2020), both including a Finnish
subsection. Opusparcus consists of candidate paraphrases automatically
extracted from the alternative translations of movie and TV show subtitles
after automatic sentence alignment. While the candidate paraphrases are
automatically extracted, a small subset of few thousand paraphrase pairs for
each language is manually annotated. TaPaCo contains candidate paraphrases
automatically extracted from the Tatoeba dataset222https://tatoeba.org/eng/,
which is a multilingual crowdsourced database of sentences and their
translations. Like Opusparcus, TaPaCo is based on language pivoting, where all
alternative translations for a same statement are collected. However, unlike
most other corpora, the candidate paraphrases are grouped into ‘sets’ instead
of pairs, and all sentences in a set are considered equivalent in meaning.
TaPaCo does not include any manual validation.

##  3 Text Selection

As discussed previously, we elect to rely on fully manual candidate extraction
as a measure against any bias introduced through heuristic candidate selection
methods. In order to obtain sufficiently many paraphrases for the person-
months spent, the text sources need to be paraphrase-rich, i.e. have a high
probability for naturally occurring paraphrases. Such text sources include for
example news headings and articles reporting on the same news, alternative
translations of the same source material, different student essays and exam
answers for the same assignment, and related questions with their replies in
discussion fora, where one can assume different writers using distinct words
to state similar meaning. For this first version of the corpus, we use two
different text sources: alternative Finnish subtitles for the same movies or
TV episodes, and headings from news articles discussing the same event in two
different Finnish news sites.

###  3.1 Alternative Subtitles

OpenSubtitles333http://www.opensubtitles.org distributes an extensive
collection of user generated subtitles for different movies and TV episodes.
These subtitles are available in multiple languages, but surprisingly often
the same movie or episode have versions in a single language, originating from
different sources. This gives an opportunity to exploit the natural variation
produced by independent translators, and by comparing two different subtitles
for a single movie or episode, there is a high likelihood of finding naturally
occurring paraphrases.

From the database dump of OpenSubtitles2018 obtained through OPUS Tiedemann
(2012), we selected all movies and TV episodes with at least two Finnish
subtitle versions. In case more versions are available, the two most lexically
differing are selected for paraphrase extraction. To filter out subtitle pairs
with low density of interesting paraphrase candidates, pairs with too high or
too low cosine similarity of TF-IDF weighted document vectors are discarded.
High similarity usually reflects identical subtitles with minor formatting
differences, while low similarity is typically caused by incorrect identifiers
in the source data. The two selected subtitle versions are then roughly
aligned using the timestamps, and divided into segments of 15 minutes. For
every movie/episode, the annotators are assigned one random such segment, the
two versions presented side-by-side in a custom tool, allowing for fast
selection of paraphrase candidates.

In total, we were able to obtain at least one pair of aligned subtitle
versions for 1,700 unique movies and TV series. While for each unique movie
only one pair of aligned subtitles is selected for annotation, TV series
comprise different episodes, dealing with the same plot and characters, and
therefore overlapping in language. After an initial annotation period, we
noticed a topic bias towards a limited number of TV series with a large number
of episodes, and decided to limit the number of annotated episodes to 10 per
each TV series in all subsequent annotation. In total, close to 3,000
different movies/episodes are used for manual paraphrase candidate extraction,
each including exactly one pair of aligned subtitles.

###  3.2 News Headings

We have downloaded news articles through open RSS feeds of different Finnish
news sites during 2017–2021, resulting in a substantial collection of news
from numerous complementary sources. For this present work, we narrow the data
down to two sources: the Finnish Broadcasting Company (YLE) and Helsingin
Sanomat (HS, English translation: Helsinki News). We align the news using a
7-day sliding window on time of publication, combined with cosine similarity
of TF-IDF-weighted document vectors induced on the article body, obtaining
article pairs likely reporting on the same event. We use the article headings
as paraphrase candidates, striving to select maximally dissimilar headings of
maximally similar articles as the most promising candidates for non-trivial
paraphrases. In practice, we used a simple grid search and human judgement to
establish the most promising region of article body and heading similarity
values.

##  4 Paraphrase Annotation

The paraphrase annotation is comprised of multiple annotation steps, including
candidate selection as described above, manual classification of candidates
based on an annotation scheme, as well as the possibility of rewriting partial
paraphrases into full paraphrases. Next, we will discuss the different
paraphrase types represented in our annotation scheme, and afterwards the
annotation workflow is discussed in a more detailed fashion.

###  4.1 Annotation Scheme

Instead of a simple yes/no (_equivalent_ or _not equivalent_) as in MRPC Dolan
and Brockett (2005) or 1–4 scale (_bad_ , _mostly bad_ , _mostly good_ and
_good_) as in Opusparcus Creutz (2018), our annotation scheme is adapted to
capture the level of paraphrasability in a more detailed fashion. Our
annotation scheme uses the base scale 1–4 similar to other paraphrase corpora,
enriched with additional subcategories (flags) for distinguishing different
types of paraphrases which would otherwise fall from the label 4 (_good_) into
label 3 (_mostly good_).

An example for each of the categories discussed below is shown in Table 2
(English translations available in Appendix A). Each candidate pair is first
evaluated in terms of the base scale numbered from 1 to 4, where:

#### Label 4

is a full (perfect) paraphrase in all reasonably imaginable contexts, meaning
one can always be replaced with the other without changing the meaning. This
ability to substitute one for the other in any context is the primary test for
_label 4_ used in the annotation.

#### Label 3

is a context dependent paraphrase, where the meaning of the two statements is
the same in the present context, but not necessarily in other contexts.

#### Label 2

is related but not a paraphrase, where there is a clear relation between the
two statements, yet they cannot be considered paraphrases.

#### Label 1

is unrelated, there being no reasonable relation between the two statements,
most likely a false positive in candidate selection.

####

If labeling a candidate pair is not possible for a reason, or giving a label
would not serve the desired purpose (e.g. wrong language or identical
statements), the example can be skipped with the _label x_.

With the base labels alone, a great number of candidate paraphrases would fail
the substitution test for label 4 and be classified _label 3_. This is
especially so for longer text segments which are less likely to express
strictly the same meaning. In order to avoid populating the _label 3_ category
with a very diverse set of paraphrases, we opt to introduce flags for finer
sub-categorization and thus support a broader range of downstream applications
of the corpus. These flags are always attached to _label 4_ (subcategories of
full paraphrases), meaning the paraphrases are not fully interchangeable due
to the specified reason, but, crucially, are context-independent, unlike
_label 3_. The possible flags are:

#### Subsumption ($>$ or $<$)

where one of the statements is more detailed and the other more general. The
relation of the pair is therefore directional, where the more detailed
statement can be replaced with the more general one in all contexts, but not
the other way around. The two common cases are one statement having additional
minor details the other omits, and one statement being ambiguous while the
other not. If there is a justification for crossing directionality (one
statement being more detailed in one aspect while the other in another
aspect), the pair falls into _label 3_ as the directional replacement test
does not hold anymore.

#### Style (s)

for tone or register difference in cases where the meaning of the two
statements is the same, but the statements differ in tone or register such
that in certain situations, they would not be interchangeable. For example, if
one statement uses pejorative language or profanities, while the other is
neutral, or one is clearly colloquial language while the other is formal. The
style flag also includes differences in the level of politeness, uncertainty,
and strength of the statements.

#### Minor deviation (i)

marks in most cases minimal differences in meaning (typically ”this” vs.
”that”) as well as easily traceable differences in grammatical number, person,
tense or such. Some applications might consider these as _label 4_ for all
practical purposes (e.g. information retrieval), while others should regard
these as _label 2_ (e.g. automatic rephrasing).

####

The flags are independent of each other and can be combined in the annotation.

Label | Statement 1 | Statement 2
---|---|---
4 | Tyrmistyttävän lapsellista! | Pöyristyttävän kypsymätöntä!
4s | Olen työskennellyt lounaan ajan. | Tein töitä koko ruokiksen.
4i | Teitäpä onnisti. | Oletpa onnekas.
4$>$ | Tein lujasti töitä niiden rahojen eteen. | Paiskin kovasti töitä.
4$<$s | Sä ruletat! Anna mennä! | Sinä olet paras, Tähkä! Anna mennä!
4is | Sä pöllit meidän kasvin! | Varastit meidän kasvit!
3 | Aion tehdä kokeen. | Aion testata sitä.
2 | Tappion kokenut Väyrynen katosi Helsingin yöhön. | Väyrynen putoamassa eduskunnasta.
| Rewrites |
Orig | Voinko palata tehtäviini? | Saanko jatkaa?
Rew | _Voinko palata tehtäviini?_ | _Saanko jatkaa tehtäviäni?_
Table 2: Example paraphrase pairs annotated with different labels and flags
(English translations available in Appendix A).

###  4.2 Annotation Workflow

Given two aligned documents as described in Section 3, an annotator first
extracts all candidate paraphrases. These can be anything between a short
phrase and several sentences long, typically being about a sentence long. The
annotators are encouraged to select as long continuous statements as possible,
nevertheless at the same time avoiding a bias towards subsumption flag by
over-extending one of the candidates. The candidate paraphrases are
subsequently transferred into a classification annotation tool. In case of
news headings, where the candidates are extracted automatically, the
candidates are introduced directly in the classification tool without any
manual extraction step.

In the classification tool, the annotator assigns a label for each candidate.
The candidate paraphrases are shown one pair at a time, and for each pair the
document context is available.

In addition to assigning a label and optional flags for a candidate pair, the
classification tool provides an option to rewrite the statements if the
classification is anything else than _label 4_ without any flags. The
annotators are instructed to rewrite the candidates in cases, where a simple
fix, for example word or phrase deletion, addition or replacement with a
synonym or changing an inflection, can be easily constructed. Rewrites must be
such that the annotated label for the rewritten example is _4_. In cases where
the rewrite would require more complicated changes or would take too much
time, the annotators are instructed to move on to the next candidate pair. One
rewrite done during the data annotation is illustrated in Table 2.

The annotators can mark unsure, difficult or otherwise interesting cases for
later discussion in daily annotation meetings. The annotators also communicate
online, for instance seeking a quick validation for a particular decision. The
work is further supported by a jointly produced 17-page annotation manual,
which is revised and extended regularly.

##  5 Corpus Statistics and Evaluation

The released corpus includes 45,663 naturally occurring paraphrases with
additional 7,909 rewrites, resulting in the total size of 53,572 paraphrase
pairs. Basic data statistics are summarized in Table 3, and label distribution
in Figure 1. Notably, the amount of candidate pairs labeled as not paraphrases
(labels 1 or 2 in our scheme) is almost non-existent, owing to the manual
candidate selection step in subtitles data from which the vast majority of the
corpus data originates. Only 5.6% of paraphrase pairs in the corpus originate
from the automated candidate selection from news data. The amount of
candidates labeled with label 1 or label x is insignificantly small, therefore
we decided to discard these from the final corpus.

Section | Examples | Rewrites | Total
---|---|---|---
Train | 36,600 | 6,239 | 42,839
Devel | 4,474 | 884 | 5,358
Test | 4,589 | 786 | 5,375
Total | 45,663 | 7,909 | 53,572
Table 3: Data sizes in our corpus. Figure 1: Labels distribution in our
corpus.

In Figure 2 we measure the density of different label combinations in the
training set conditioned on cosine similarity of paraphrase pairs based on TF-
IDF weighted character n-grams of lengths 2–4. Up to cosine similarity of 0.5
the most common labels are evenly represented, while the prevalence of label 4
increases throughout the range and dominates the sparsely populated range of
similarities over 0.8.

Figure 2: Density of different labels in the training set conditioned on
cosine similarity of the paraphrase pairs.

###  5.1 Annotation Quality

After the initial annotator training phase most of the annotation work is
carried out as single annotation. In order to monitor annotation consistency,
double annotation batches are assigned regularly. In double annotation, one
annotator first extracts the candidate paraphrases from the aligned documents,
but later on these candidates are assigned to two different annotators, who
annotate the labels for these independently from each other. Next, these two
individual annotations are merged and conflicting labels are resolved together
with the whole annotation team in a meeting. These consensus annotations
constitute a gold standard against which individual annotators can be
measured.

A total of 1,175 examples are double annotated (2.5% of the data444During the
initial annotator training double annotation was used extensively; this
annotator training data is not included in the released corpus.). Most of
these are annotated by exactly two annotators, however, some examples may
include annotations from more than two annotators, and thus the total amount
of individual annotations for which the gold standard label exists is 2,513.
We measure the agreement of individually annotated examples against the gold
standard annotations in terms of accuracy, i.e. the proportion of individually
annotated examples with correctly assigned label.

The overall accuracy is 68.7% when the base label (_labels 1–4_) as well as
all additional flags are taken into consideration. When discarding the least
common flags _s_ and _i_ and evaluating only base labels and directional
subsumption flags, the overall accuracy is 72.9%. To compare the observed
agreement to previous studies on paraphrase annotation, the Opusparcus
annotation agreement is approximately 64% on Finnish development set and 67%
on test set (calculated from numbers in Table 4 and Table 5 in Creutz (2018)).
The Opusparcus uses an annotation scheme with four labels, similar to our base
label scheme. In MRPC, the reported agreement score is 84% on a binary
paraphrase-or-not scheme. While direct comparison is difficult due to the
different annotation schemes and label distributions, the figures show that
the observed agreement seem to be roughly within the same range with agreement
numbers seen in related works.

In addition to agreement accuracy, we calculate two versions of Cohen’s kappa,
a metric for inter-annotator agreement taking into account the possibility of
agreement occurring by chance. First we measure kappa agreement of all
individual annotations against the gold standard, an approach typical in
paraphrase literature. This kappa is 0.62, indicating substantial agreement.
Additionally, we measure the Cohen’s kappa between each pair of annotators.
The weighted average kappa over all annotator pairs is 0.41 indicating
moderate agreement. Both are measured on full labels. When evaluating only on
base labels and directional subsumption flags, these kappa scores are 0.65 and
0.45, respectively.

###  5.2 Corpus Comparison

Figure 3: Comparison of paraphrase length distributions in terms of tokens per
paraphrase.

Figure 4: Comparison of paraphrase pair cosine similarity distributions.

We compare the distribution of paraphrase lengths and lexical similarity with
the two Finnish paraphrase candidate corpora, Opusparcus and TaPaCo, as the
reference. Direct comparison is complicated by several factors. Firstly, both
Opusparcus and TaPaCo consist primarily of automatically extracted paraphrase
candidates, Opusparcus having only small manually curated development and test
sections, and TaPaCo being fully uncurated. Secondly, the small manually
annotated sections of Opusparcus are sampled to emphasize lexically dissimilar
pairs, and therefore not representative of the characteristics of the whole
corpus. We therefore compare with the fully automatically extracted sections
of both Opusparcus and TaPaCo. For our corpus, we discard the small proportion
of examples of base labels 1 and 2, i.e. not paraphrases. Another important
factor to consider is that the proportion of false candidates in the
automatically extracted sections of Opusparcus and TaPaCo is unknown, further
decreasing comparability: the characteristics of false and true candidates may
differ substantially, false candidates for example likely being on average
more dissimilar in terms of lexical overlap than true candidates.

For each corpus, we sample 12,000 paraphrase pairs. For our corpus, we
selected a random sample of true paraphrases (_label 3_ or higher) from the
train section. For TaPaCo, the sample covers all paraphrase candidates from
the corpus, however with the restriction of taking only one, random pair from
each ‘set’ of paraphrases. For Opusparcus, which is sorted by a confidence
score in descending order, the sample was selected to contain the most
confident 12K paraphrase candidates.555When we repeated the length analysis
with a sample of 480K most confident pairs, the length distribution and
average length remained largely unchanged, while the similarity distribution
became close to flat. Without manual annotation, it is hard to tell the reason
for this behavior.

In Figure 3 the length distribution of paraphrases in terms of tokens is
measured for the abovementioned samples. Although the majority of paraphrases
are rather short in all three corpora, we see that our corpus includes
considerably higher proportion of longer paraphrases. The average number of
tokens in our corpus is 8.3 tokens per paraphrase, while it is 5.6 in TaPaCo
and 3.6 in Opusparcus candidates.

In Figure 4 the paraphrase pair cosine similarity distribution is measured
using TF-IDF weighted character n-grams of length 2–4. While both TaPaCo and
Opusparcus lean towards higher similarity candidates, the distribution of our
corpus is more balanced including a considerably higher proportion of pairs
with low lexical similarity.

##  6 Paraphrase Classification Baseline

In order to establish a baseline classification performance on the new
dataset, we train a classifier based on the FinBERT model Virtanen et al.
(2019). Each paraphrase pair of statements A and B is encoded as the sequence
[CLS] A [SEP] B [SEP], where [CLS] and [SEP] are the special marker tokens of
the BERT model. Subsequently, the output embeddings of the three special
tokens are concatenated together with the averaged embeddings of the tokens in
A and B. These five concatenated embeddings are then propagated into four
decision layers: one for the base label _2/3/4_ , one for the subsumption flag
_$ <$/$>$/none_, and one for each the binary flag _s_ and _i_. Since the flags
only apply to base label _4_ , no gradients are applied to these layers for
examples with base labels _2_ and _3_. We have explored also other BERT-based
architectures, such as basing the classification on the [CLS] embedding only
as is customary, and having a single classification layer comprising all
possible base label and flag combinations. These resulted in a consistent drop
in prediction accuracy, and we did not pursue them any further.

The baseline results are listed in Table 4 showing that per-class F-score
ranges between 38–71%, strongly correlated with the number of examples
available for each class. When interpreting the task as a pure multi-class
classification, i.e. when counting all possible combinations of base label and
flags as their own class, the accuracy is 54% with majority baseline being
34.3%, and the annotators’ accuracy 68.7%. The model thus positions roughly to
the mid-point between the trivial majority baseline, and human performance.

Label | Prec | Rec | F-score | Support
---|---|---|---|---
2 | 50.9 | 31.2 | 38.7 | 93
3 | 57.7 | 31.9 | 41.1 | 990
4 | 66.2 | 78.2 | 71.7 | 2149
4< | 52.8 | 53.5 | 53.2 | 1007
4> | 52.6 | 56.1 | 54.3 | 1136
i | 51.5 | 36.5 | 42.7 | 329
s | 51.4 | 28.9 | 37.0 | 249
W. avg | 52.9 | 54.0 | 52.2 |
Acc |  |  | 54.0 |
Table 4: Classification performance on the test set, when the base label and
the flags are predicted separately. In the upper section, we merge the
subsumption flags with the base class prediction, but leave the _i_ and _s_
separated. The rows _W. avg_ and _Acc_ on the other hand refer to performance
on the complete labels, comprising all allowed combinations of base label and
flags. _W. avg_ is the average of P/R/F values across the classes, weighted by
class support. _Acc_ is the accuracy.

##  7 Discussion and Future Work

In this work, we set out to build a paraphrase corpus for Finnish that would
be (a) in the size category allowing deep model fine-tuning and (b) manually
gathered maximizing the chance of finding more non-trivial, longer paraphrases
than would be possible with the traditional automatic candidate extraction.
The annotation so far took 14 person-months and resulted in little over 50,000
manually classified paraphrases. We have demonstrated that, indeed, the corpus
has longer, more lexically dissimilar paraphrases. Building such a corpus is
therefore shown feasible and presently it is likely the largest manually
annotated paraphrase dataset for any language, naturally at the inevitably
higher data collection cost. The manual selection is only feasible for texts
rich in paraphrase, and the domains and genres covered by the corpus is
necessarily restricted by this condition.

In our future work, we intend to extend the manually annotated corpus, ideally
roughly double its present size. We expect the pursued data size will allow us
to build sufficiently accurate models, both in terms of embedding and pair
classification, to gather further candidates automatically at a level of
accuracy sufficient to support down-stream applications. We are also
investigating further text sources, especially parallel translations outside
of the present subtitle domain. The additional flags in our annotation scheme,
as well as the nearly 10,000 rewrites allow for interesting further
investigations in their own right.

##  8 Conclusions

In this paper we presented the first entirely manually annotated paraphrase
corpus for Finnish including 45,663 naturally occurring paraphrases gathered
from alternative movie or TV episode subtitles and news headings. Further
7,909 hand-made rewrites are provided, turning context-dependent paraphrases
into perfect paraphrases whenever possible. The total size of the released
corpus is 53,572 paraphrase pairs of which 98% are manually classified to be
at least paraphrases in their given context if not in all contexts.

Additionally, we evaluated the advantages and costs of manual paraphrase
candidate selection from two ‘parallel’ but monolingual documents. We
demonstrated the approach on alternative subtitles showing the technique being
feasible for high quality candidate selection yielding sufficient amount of
paraphrase candidates for the given annotation effort. We have shown the
candidates to be notably longer and less lexically overlapping than what
automated candidate selection permits.

The corpus is available at github.com/TurkuNLP/Turku-paraphrase-corpus under
the CC-BY-SA license.

## Acknowledgments

We gratefully acknowledge the support of European Language Grid which funded
the annotation work. Computational resources were provided by _CSC — the
Finnish IT Center for Science_ and the research was supported by the Academy
of Finland. We also thank Sampo Pyysalo for fruitful discussions and feedback
throughout the project.

## References

  * Altheneyan and Menai (2019) Alaa Saleh Altheneyan and Mohamed El Bachir Menai. 2019.  Evaluation of state-of-the-art paraphrase identification and its application to automatic plagiarism detection.  _International Journal of Pattern Recognition and Artificial Intelligence_ , 34.
  * Bhagat and Hovy (2013) Rahul Bhagat and Eduard Hovy. 2013.  Squibs: What is a paraphrase?  _Computational Linguistics_ , 39(3):463–472.
  * Creutz (2018) Mathias Creutz. 2018.  Open Subtitles paraphrase corpus for six languages.  In _Proceedings of the Eleventh International Conference on Language Resources and Evaluation (LREC 2018)_.
  * Dolan and Brockett (2005) William B. Dolan and Chris Brockett. 2005.  Automatically constructing a corpus of sentential paraphrases.  In _Proceedings of the Third International Workshop on Paraphrasing (IWP 2005)_.
  * Dong et al. (2021) Qingxiu Dong, Xiaojun Wan, and Yue Cao. 2021.  ParaSCI: A large scientific paraphrase dataset for longer paraphrase generation.  _arXiv preprint arXiv:2101.08382_.
  * He et al. (2020) Yun He, Zhuoer Wang, Yin Zhang, Ruihong Huang, and James Caverlee. 2020.  PARADE: A new dataset for paraphrase identification requiring computer science domain knowledge.  In _Proceedings of the 2020 Conference on Empirical Methods in Natural Language Processing (EMNLP)_ , pages 7572–7582.
  * Lan et al. (2017) Wuwei Lan, Siyu Qiu, Hua He, and Wei Xu. 2017.  A continuously growing dataset of sentential paraphrases.  In _Proceedings of the 2017 Conference on Empirical Methods in Natural Language Processing (EMNLP)_ , pages 1224–1234.
  * Madnani and Dorr (2010) Nitin Madnani and Bonnie J. Dorr. 2010.  Generating phrasal and sentential paraphrases: A survey of data-driven methods.  _Computational Linguistics_ , 36(3):341–387.
  * Mehdizadeh Seraj et al. (2015) Ramtin Mehdizadeh Seraj, Maryam Siahbani, and Anoop Sarkar. 2015.  Improving statistical machine translation with a multilingual paraphrase database.  In _Proceedings of the 2015 Conference on Empirical Methods in Natural Language Processing (EMNLP)_ , pages 1379–1390.
  * Scherrer (2020) Yves Scherrer. 2020.  TaPaCo: A corpus of sentential paraphrases for 73 languages.  In _Proceedings of The 12th Language Resources and Evaluation Conference (LREC 2020)_ , pages 6868–6873.
  * Soni and Roberts (2019) Sarvesh Soni and Kirk Roberts. 2019.  A paraphrase generation system for EHR question answering.  In _Proceedings of the 18th BioNLP Workshop and Shared Task_ , pages 20–29.
  * Tiedemann (2012) Jörg Tiedemann. 2012.  Parallel data, tools and interfaces in OPUS.  In _Proceedings of the Eighth International Conference on Language Resources and Evaluation (LREC 2012)_ , pages 2214–2218.
  * Virtanen et al. (2019) Antti Virtanen, Jenna Kanerva, Rami Ilo, Jouni Luoma, Juhani Luotolahti, Tapio Salakoski, Filip Ginter, and Sampo Pyysalo. 2019.  Multilingual is not enough: BERT for Finnish.  _arXiv preprint arXiv:1912.07076_.

##  Appendix A English Translation of Table 2

Label | Statement 1 | Statement 2
---|---|---
4 | Shockingly childish! | Astoundingly immature!
4s | I have worked for the duration of lunch. | I worked through the whole chowtime.
4i | You guys got lucky, didn’t you. | Aren’t you fortunate.
4$>$ | I worked so hard for the money. | I put so much effort into work.
4$<$s | You rule! Come on, dude! | You are the best, Tähkä! Come on!
4is | You nicked our plant! | You stole our plants!
3 | I intend to make an experiment. | I am going to test it.
2 | Defeated Väyrynen vanished into the | Väyrynen is losing his seat in the parliament.
| Helsinki night |
| Rewrites |
Orig | Can I get back to my assignments? | Can I continue?
Rew | _Can I get back to my assignments?_ | _Can I continue working on my assignments?_
Table 5: English translations for annotation examples in Table 2.