E_mail<EMAIL_ADDRESS>
Dipartimento di Ingegneria e Scienze dell’Informazione e Matematica,
Università degli Studi dell’Aquila, via Vetoio, 67100 L’Aquila, Italy.

E_mail<EMAIL_ADDRESS>
Dipartimento di Ingegneria e Scienze dell’Informazione e Matematica,
Università degli Studi dell’Aquila, via Vetoio, 67100 L’Aquila, Italy.

E_mail<EMAIL_ADDRESS>
Dipartimento di Scienze Matematiche, Politecnico di Torino, Corso Duca degli
Abruzzi 24, 10129 Torino, Italy
INFN, Sezione di Torino, Via Pietro Giuria 1, 10125 Torino, Italy
ORCID: 0000-0002-4223-6279

# Finite reservoirs and irreversibility corrections to Hamiltonian systems
statistics

Matteo Colangeli  Antonio Di Francesco  Lamberto Rondoni

###### Abstract

We consider several Hamiltonian systems perturbed by external agents, that
preserve their Hamiltonian structure. We investigate the corrections to the
canonical statistics resulting from coupling such systems with possibly large
but finite reservoirs, and from the onset of processes breaking the time
reversal symmetry. We analyze exactly solvable oscillators systems, and
perform simulations of relatively more complex ones. This indicates that the
standard statistical mechanical formalism needs to be adjusted, in the ever
more investigated nano-scale science and technology. In particular, the
hypothesis that heat reservoirs be considered infinite and be described by the
classical ensembles is found to be critical when exponential quantities are
considered, since the large size limit may not coincide with the infinite size
canonical result. Furthermore, process-dependent emergent irreversibility
affects ensemble averages, effectively frustrating, on a statistical level,
the time reversal invariance of Hamiltonian dynamics, that is used to obtain
numerous results.

Keywords: Jarzynski equality, nonequilibrium process, fluctuation relations,
finite size effects

##  1 Introduction

The validity of the canonical ensemble is universally accepted to compute
macroscopic quantities of systems in equilibrium at a given temperature $T$,
as averages of phase space functions. The corresponding formalism assumes that
heat reservoirs are infinitely large, and that measurement times are
exceedingly longer than the characteristic times of the microscopic events.
Such mathematical idealizations yield a highly successful theory describing a
vast range of macroscopic phenomena. The separation between microscopic and
macroscopic scales is indeed sufficiently wide for calculations of quantities
of thermodynamic interest. Nevertheless, there are various reasons for
investigating the applicability of the canonical framework to non-standard
observables. For instance, exponentials of microscopically expressed variables
appear in Bennett’s formulae for the free energy [1], in Widom’s relation [2],
Zwanzig’s relation [3], and in the more recent Jarzynski [4] and Crooks
relations [5]. Furthermore, current science and technology deal with small
systems and fast processes, as well as with quantities not immediately
interpretable in thermodynamic terms, as in the case of anomalous enegy
transport [6, 7, 8]. Therefore, finite size effects and lack of ergodicity may
turn important. Indeed, standard thermodynamic properties of macroscopic
objects only require a proper characterization of the bulk of the relevant
probability distributions, not of their tails. On the other hand, an accurate
characterization of the tails of the relevant probability distributions
becomes necessary when dealing with observables that get a substantial
contribution from such tails. Then, the fact that thermal baths are
necessarily finite and that experiments may last very short times may require
particular attention.

In this work we take the quantity used in the Jarzynski Equality (JE) as a
paradigmatic example of topical non-standard observables. It is worth
recalling that the time reversal symmetry of the microscopic dynamics [9, 10,
11, 12] is essential for the derivation of the JE, which belongs to a class of
results, known as Fluctuation Relations, strongly relying on the time
reversibility of the microscopic dynamics, see e.g. [13, 14, 15]. More
generally, the time reversal symmetry turns out being a standard ingredient of
a large variety of statistical mechanical results, including the Onsager
Reciprocal Relations [16, 17], Fluctuation-Dissipation Theorem and the Green-
Kubo relations [18, 19, 20], and applications to magnetic systems [21]. In
works such as Ref.[22] it was found that certain nanoscopic Hamiltonian
systems violate the JE, although formally amenable to analysis within the
canonical framework, that yields the JE as an exact relation [4]. In the case
of Ref.[22], the failure was caused by the emergence of irreversibility due to
a process dependent nonequilibrium effect, and not to the large number of
degrees of freedom. At the same time, highly nonequilibrium processes do not
prevent the validity of the JE in e.g. 1-dimensional systems described by an
overdamped Langevin equation [23]. That this may be the case is clear in the
words of, e.g Fermi [24] or Callen [25], who state, in practice, that the
ensembles work if the observation times suffice for the observables of
interest to have thouroughly explored their range. Khinchin then adds that
this is easy to obtain, for the observables of interest typically have a small
range [26]. In all instances, the state of the system is required to be
stationary, or very slowly evolving with respect to the observation times.

The above considerations are topical, given the rapid development of bio- and
nano-technologies, which deal with small systems. Apart from being small, such
systems are often briskly driven by external agents, so that thermal baths
(even if effectively infinite) only express limited energy, and the
deterministic thermodynamic laws must often be replaced by statistical laws.
Certainly, some experiments of bio- and nano-technological interest
intentionally take very short times, so that only a small part of a thermal
bath is effectively involved. This poses the question, when computing ensemble
averages, about proper approaches to the finiteness of the bath or, in other
terms, to the restriction to finite subsets of phase space.

In this paper we thus analyze the finite size effects on the statistics
concerning simple Hamiltonian systems, subjected to various external drivings.
We start by briefly reviewing the derivation of the JE in Sec. 2. In Sec. 3,
three simple mechanical models are introduced to illustrate the onset of
finite size effects that lead to violations of the JE, highlighting some of
the limitations of the canonical statistic. In particular it is shown that the
speed of of the protocol or the frequency of periodic drivings resonating with
the system proper frequency may wildly enhance the protocol dependence,
violating of up to 110% the JE. We also highlight the fact that analagous
results are obtained for infinite baths at small temperatures.

In Sec. 4, we consider a model mimicking the adiabatic expansion of an ideal
gas, and also describe the validity of the JE in the presence of a protocol-
dependent device concluding that the occurrence of an irreversible phenomenon
(such as the free expansion of a gas) can invalidate the statistical
description of a particle system through the canonical formalism.

Conclusions are drawn in Sec. 5, where we also anticipate future developments.
The Appendices give the details of some analytical calculations reported in
the main text.

##  2 Derivation of the Jarzynski equality

A well-known example involving both the canonical ensemble and exponential
variables is the Jarzynski Equality (JE), which offers a useful playground to
highlight the role of finite size effects on the statistics of thermodynamics
quantities in the canonical framework. The JE has been derived for both
stochastic, and deterministic systems. We focus on the second, which concerns
a system S made of $N$ particles, initially in equilibrium with a bath B at
temperature $T$. The system may interact with an environment, E, also
initially in equilibrium with B. The Hamiltonian of system and environment,
denoted by S+E, is assumed to take the following form:

${\cal
H}(\mathbf{x},\mathbf{v};\lambda)=H_{S}(x_{S},v_{S};\lambda)+H_{E}(x_{E},v_{E})+h_{I}(\mathbf{x},\mathbf{v})$
(2.1)

where $\lambda$ is a parameter controlled by an external agent,
$(\mathbf{x},\mathbf{v})=(x_{S},v_{S},x_{E},v_{E})$ are the position and
velocities of S and E, as indicated by the subscripts, $H_{S}$, $H_{E}$ and
$h_{I}$ are, respectively, the energy of S, the energy of E and the energy of
their interaction. The initial distributions of coordinates and momenta of
S+E, which is in equilibrium at temnperature $T$, is given by the canonical
ensemble:

$P_{0}(\Gamma)={1\over Z_{0}}e^{-\beta{\cal
H}(\Gamma;A)}\,,\quad\beta=\cfrac{1}{k_{{}_{B}}T}$ (2.2)

where $k_{B}$ is the Boltzmann’s constant, $\Gamma=(\mathbf{x},\mathbf{v})$ is
one configuration of S+E, and $Z_{0}$ is the initial canonical partition
function.

At time $t=0$, this system is isolated from the bath, and driven by an
external agent that modifies the parameter $\lambda$. This is done many times,
repeating the same protocol $\lambda:[0,\tau]\rightarrow\mathbb{R}$ over a
given finite time $\tau$, changing the parameter from its initial value
$\lambda(0)=A$, to its final value $\lambda(\tau)=B$. Each time, a different
initial condition is taken at random, according to the canonical distribution
(1.43), and the following quantity, called work, is computed [4]:

$W_{J}(\Gamma_{0})=\int_{0}^{\tau}\,{\partial
H\over\partial\lambda}\,\dot{\lambda}\,\mbox{d}t={\cal
H}(\Gamma_{\tau}(\Gamma_{0});B)-{\cal H}(\Gamma_{0};A)$ (2.3)

where $\Gamma_{\tau}(\Gamma_{0})$ is the phase reached in the time $\tau$
starting from the initial condition $\Gamma_{0}$. Because the protocol
$\lambda(t)$ is fixed, the dynamics are deterministic, and the value of the
work depends only on the initial condition. However, the initial conditions
change randomly, yielding a different value of $W_{J}$ for each realization of
the process, and effectively making it a random variable. In this setting, the
following relation, known as Jarzynski equality, was obtained: [4]:

$\Big{\langle}e^{-\beta W_{J}}\Big{\rangle}_{0}=e^{-\beta\Delta F}$ (2.4)

where $\langle\cdot\rangle_{0}$ is the canonical ensemble average obtained
from $P_{0}$, and $\Delta F=F_{B}-F_{A}$ is the equilibrium free energy
difference between the equilibrium canonical state with parameter $\lambda=B$
and the one with parameter $\lambda=A$, both at temperature $T$. One of the
most striking aspects of the JE, that is a direct effect of the canonical
ensemble, is that it does not depend on the protocol. This sounds at odds with
the fact that physical theories have a range of applicability limited by space
and time constraints, outside of which a different description must be
adopted. On the other hand, it depends on the validity of the canonical
ensemble whose applicability boundaries are not known, in general, especially
if involving non standard quantities. Understanding the role of the canonical
ensemble is important in general, not just in relation to the JE. We will see
that the quantity in the right hand side of Eq.(2.4), depends on the protocol,
if the ensemble does not extend to infinity. Note that the form of the
probability distribution properly describing the effect of finite environments
are not known in general, but the finite size effects can be evidenced on any
distribution. In the concluding remarks we address this issue.

##  3 Models and methods

Below, we investigate possible finite size effects for several different
systems. In particular, we analyze three simple harmonic oscillators models,
perturbed from their equilibrium states. The perturbation is applied by
harmonic springs, whose center of force moves according to deterministic rules
$\lambda(t)$. The first model consists of a single oscillator, playing the
role of S, with $\lambda(t)=\ell t$, $t\in[0,\tau]$, where $\ell$ and $\tau$
are constants that can be varied, in such a way that the initial and final
values of $\lambda$ do not change: $\lambda(0)=A$ and $\lambda(\tau)=B$. In
the second model, the protocol is changed to $\lambda(t)=\sin(\gamma t)$.
Being periodic in time, this protocol yields different phenomena when the
frequency $\gamma$ is changed, like resonances that affect $W_{J}$ and,
consequently the JE. Both, the first and the second case do not have any
environment E or, equivalently, the interaction energy vanishes: $h_{I}=0$.
The third model we consider has two oscillators, one of which is taken to be
the system S and the other the environment E. As the theory requires, only S
is subjected to a time dependent perturbation.

###  3.1 Single oscillator under linear protocol

Take a 1D system made of a single harmonic oscillator with rest position in
$x=0$, that is driven by a moving harmonic trap, centered in $\lambda(t)=\ell
t$, where $\ell$ is a positive constant, and $t\in[0,\tau]$. The initial value
of $\lambda$ is given by $A=\lambda(0)=0$, and let its final value be denoted
by $B=\lambda(\tau)=\ell\tau$, with $B$ fixed. To explore the effect of
modifying the speed of the protocol, we vary $\ell$ and $\tau$, so that $B$ is
fixed. Let the oscillator mass be $m$, and its momentum $p=mv$, where $v$ is
the velocity. Then, the motion is determined by the following time dependent
Hamiltonian:

$\displaystyle{\cal H}(x,v;t)$ $\displaystyle=$ $\displaystyle{p^{2}\over
2m}+{k_{p}\over 2}x^{2}+{k_{D}\over 2}\left(\lambda-x\right)^{2}={p^{2}\over
2m}+{k\over 2}x^{2}+{k_{D}\ell\over 2}\left(\ell t^{2}-2xt\right)~{},$ (3.5)

where $k_{p}$ is the elastic constant of the spring with rest position in
$x=0$, $k_{D}$ the elastic constant of the moving trap, and $k=k_{D}+k_{p}$.
The equation of motion consequently takes the form:

$\ddot{x}=-\omega^{2}x+{k_{D}\over m}\ell t~{},\quad\mbox{with
}~{}x(0)=x_{0}\,,~{}v(0)=v_{0}$ (3.6)

where we introduced the natural frequency of the oscillator
$\omega=\sqrt{k/m}$. In this case, the work $W_{J}$ is expressed by:

$\displaystyle W_{J}=\int_{0}^{\tau}\,k_{D}\left(\ell
t-x\right)\ell\,\mbox{d}t={k_{D}\ell^{2}\tau^{2}\over
2}-k_{D}\ell\int_{0}^{\tau}\,x(t;x_{0},v_{0})\,\mbox{d}t$ $\displaystyle\hskip
20.0pt=k_{D}B\left[{B\over
2}-{1\over\tau}\int_{0}^{\tau}\,x(t;x_{0},v_{0})\,\mbox{d}t\right]$ (3.7)

where the oscillator position is expressed by:

$x(t;x_{0},v_{0})=x_{0}\cos\omega t+{v_{0}-\ell k_{D}/k\over\omega}\sin\omega
t+{\ell k_{D}\over k}t$ (3.8)

Then, performing the integration in expression (3.7), one obtains:

$W_{J}(\ell;x_{0},v_{0})=k_{D}B\left[{B\over 2}\left(1-{k_{D}\over
k}\right)-{x_{0}\ell\over B\omega}\sin\omega{B\over\ell}+\left({p_{0}\ell\over
Bk}-{\ell^{2}k_{D}\over
Bk\omega^{2}}\right)\left(\cos\omega{B\over\ell}-1\right)\right]~{},$ (3.9)

where $B$ is fixed, while the protocol speed $\ell$ can be varied. Although
$\exp(-\beta W_{J})$ depends on $\ell$, its average with respect to the
initial canonical ensemble, $P_{0}$, does not. Given $\ell\in(0,\infty)$, one
has:

$\Big{\langle}e^{-\beta
W_{J,\ell}}\Big{\rangle}_{0}=\exp\left\\{-\beta{k_{D}k_{p}B^{2}\over
2k}\right\\}\;,$ (3.10)

which does not depend on the speed of the protocol, as the Jarzynski theory
predicts. Explicit calculations are reported in the Appendix A.

In the case in which the environment is bounded and the bath can only express
a finite energy, the corresponding probability density is truncated at a given
distance $L$ from the rest position of the oscillator, and at a maximum
momentum $M$. For the sake of argument, we assume that the form of the finite
support distribution is the canonical one, truncated and normalized, and that
the two bounds $L$ and $M$ do not depend on each other. After all, the
classical ensembles constitute a most successful postulate of statistical
mechanics that, however, only seldom can be derived from the particles
dynamics. Morevoer the resulting distributions are truncated Gaussians, hence
mitigte the effects of truncation. Then, suppose we have:

$P_{0}(x,p)=\cfrac{1}{Z_{0}(L,M)}\left\\{\begin{array}[]{lcr}e^{-\beta{(kx^{2}+p^{2}/m)/2}}&\mbox{if}&|x|\leq
L~{}\mbox{ and }~{}|p|\leq M\\\\[7.0pt] 0&\mbox{if}&|x|>L~{}\mbox{ or
}~{}|p|>M\end{array}\right.$ (3.11)

with $Z_{0}(L,M)$ a normalizing factor. In this case, one obtains:

$\Big{\langle}e^{-\beta W_{J,\ell}}\Big{\rangle}_{0;L,M}=I_{exp}\cdot
I_{x}\cdot I_{p}$ (3.12)

where $I_{exp}$ represents the infinite size result, that does not depend on
$\ell$, while the finite size correction factors $I_{x}$ and $I_{p}$ do depend
on $\ell$, hence on the protocol. The explicit expressions of
$I_{exp},I_{x},I_{p}$, along with the detailed calculations leading to Eq.
(3.12), are deferred to the Appendix A. This result shows that for fixed
$\ell$ and $m$, sufficiently large $L$ and $M$ exist such that the infinite
size result is recovered; indeed $I_{x}$ and $I_{p}$ both tend to $1$, if
$L,M$ grow at fixed $\ell$. However, for fixed $L$ and $M$, sufficiently large
$\ell$, i.e. a sufficiently fast protocol, together with a large enough value
of the product $\omega B$, or sufficiently large $m$, yield $I_{x},I_{p}<1$,
i.e. $\exp(-\beta W_{J})\rangle_{0;L,M}<\langle\exp(-\beta W_{J})\rangle_{0}$.
The term $I_{p}$ is particularly sensitive to variations of $m$, because the
argument of the error function on the right of its numerator may even turn
negative, if $m$ is sufficiently large. In any event, the left hand side of
the JE is protocol dependent, if the ensemble is finitely supported. While at
odds with the infinite bath result, this is in accord with the fact that too
fast protocols (e.g. comparable with microscopic rates) require a specifically
developed approach.

###  3.2 Single Oscillator with Periodic Forcing

Let now the single oscillator be driven by a moving harmonic trap centered in
$\lambda(t)=\sin{\gamma t}$, where $\gamma={2\pi}/{T}$, and $T$ is the period
of the center of force of the trap. Take $t\in[0,\tau]$, $A=\lambda(0)=0$, and
$B=\lambda(\tau)=\sin{\gamma\tau}$. If the final value of $\lambda$ is fixed,
as in the previous subsection, different $\gamma$ correspond to faster or
slower protocols, that last a time $\tau=\arcsin({B})/\gamma$. The time
dependent Hamiltonian now takes the form

$\displaystyle{\cal H}(x,v;t)$ $\displaystyle=$ $\displaystyle{p^{2}\over
2m}+{k_{p}\over 2}x^{2}+{k_{D}\over 2}\left(\lambda-x\right)^{2}$ (3.13)
$\displaystyle=$ $\displaystyle{p^{2}\over 2m}+{k_{p}\over 2}x^{2}+{k_{D}\over
2}\left(\sin{\gamma t}-x\right)^{2}={p^{2}\over 2m}+{k\over
2}x^{2}+{k_{D}\over 2}\left(\sin^{2}\gamma t-2x\sin{\gamma t}\right)~{},$
(3.14)

where $k=k_{D}+k_{p}$. In this case, the Jarzynski work $W_{J}$ is given by:

$\displaystyle W_{J}=\int_{0}^{\tau}\,{\partial
H\over\partial\lambda}\,\dot{\lambda}\,\mbox{d}t=$ (3.15) $\displaystyle\hskip
20.0pt={k_{D}\over 4}\left(1-\cos
2\gamma\tau\right)-k_{D}\gamma\int_{0}^{\tau}x(t;x_{0},v_{0})\cos\left(\gamma
t\right)\,\mbox{d}t$ (3.16)

Given the Hamiltonian (3.14), the equation of motion for this system is:

$\displaystyle\ddot{x}=-\omega^{2}x+{k_{D}\over m}\sin\gamma
t~{},\quad\mbox{with i.c.\ }~{}x(0)=x_{0}\,,~{}v(0)=v_{0}$ (3.17)

where $\omega=\sqrt{k/m}$. For $\gamma\neq\omega$, one obtains:

$\displaystyle x(t;x_{0},v_{0})={k_{D}/m\over\omega^{2}-\gamma^{2}}\sin\gamma
t+{1\over\omega}\left(v_{0}-{\gamma
k_{D}/m\over\omega^{2}-\gamma^{2}}\right)\sin\omega t+x_{0}\cos\omega t$
(3.18)

and the work takes the form:

$\displaystyle\hskip-60.0ptW_{J}(\gamma;x_{0},v_{0})={k_{D}\over
4}\left(1-\cos 2\gamma\tau\right)-k_{D}\gamma\int_{0}^{\tau}\cos\left(\gamma
t\right)\times$ (3.19) $\displaystyle\hskip
30.0pt\left[{k_{D}/m\over\omega^{2}-\gamma^{2}}\sin\gamma
t+{1\over\omega}\left(v_{0}-{\gamma
k_{D}/m\over\omega^{2}-\gamma^{2}}\right)\sin\omega t+x_{0}\cos\omega
t\right]\,\mbox{d}t$ (3.20)

Solving the integral on the right, we finally get:

$\displaystyle\hskip-50.0ptW_{J}(\gamma;x_{0},v_{0})={k_{D}\over
4}\left(1-{k_{D}/m\over\omega^{2}-\gamma^{2}}\right)\left(1-\cos
2\gamma\tau\right)+$ (3.21) $\displaystyle\hskip
20.0pt-{k_{D}\gamma\over\omega^{2}-\gamma^{2}}\left(v_{0}-{k_{D}\gamma/m\over\omega^{2}-\gamma^{2}}\right)\left(1-{\gamma\over\omega}\sin\gamma\tau\sin\omega\tau-\cos\gamma\tau\cos\omega\tau\right)+$
(3.22) $\displaystyle\hskip
20.0pt+\,x_{0}{k_{D}\gamma\omega\over\omega^{2}-\gamma^{2}}\left({\gamma\over\omega}\sin\gamma\tau\cos\omega\tau-\cos\gamma\tau\sin\omega\tau\right)$
(3.23)

This quantity can now be multiplied by $-\beta$, exponentiated and averaged
over all the initial conditions $(x_{0},v_{0})$. In the case of the full
canonical ensemble, one obtains a result that does not depend on $\gamma$,
when $A$, $\tau$ and consequently $B$ are fixed. If, on the other hand, the
probability density is expressed by Eq.(3.11), one finds:

$\Big{\langle}e^{-\beta W_{J}}\Big{\rangle}_{0;L,M}=I_{exp}\cdot I_{x}\cdot
I_{p}\;,$ (3.24)

where the explicit expressions of $I_{exp},I_{x},I_{p}$ are given in the
Appendix B. The resonance, corresponding to $\gamma=\omega$, must be treated
separately, since the solution of the equation of motion (3.17) takes the
form:

$\displaystyle x(t;x_{0},v_{0})=x_{0}\cos\omega t+{v_{0}\over\omega}\sin\omega
t+{k_{D}/m\over 2\omega^{2}}\left(\sin\omega t-\omega t\cos\omega t\right)$
(3.25)

The Jarzynski work is now expressed by:

$\displaystyle\hskip-40.0ptW_{J}(\tau;x_{0},v_{0})={k_{D}^{2}/m\over
8}\tau^{2}+{k_{D}^{2}/m\over 8\omega}\tau\sin 2\omega\tau+{k_{D}\over
4}\left(1-{3\over 4}{k_{D}/m\over\omega^{2}}\right)\left[1-\cos
2\omega\tau\right]$ $\displaystyle\hskip 60.0pt-x_{0}{k_{D}\over
2}\left[\omega\tau+{1\over 2}\sin 2\omega\tau\right]-v_{0}{k_{D}\over
4\omega}\left[1-\cos 2\omega\tau\right]$ (3.26)

and its finite energy ensemble average can again be written as:

$\Big{\langle}e^{-\beta
W_{J}}\Big{\rangle}_{0;L,M}^{\rm(res)}=I_{exp}^{\rm(res)}\cdot
I_{x}^{\rm(res)}\cdot I_{p}^{\rm(res)}$ (3.27)

where:

$\displaystyle I_{exp}^{\rm(res)}=\exp\left\\{{\beta k_{D}\over
2}\left({k_{D}\over k}-1\right)\sin^{2}\omega\tau\right\\}$ (3.28)

$\displaystyle I_{x}^{\rm(res)}={1\over 2~{}\mbox{erf}\left(\sqrt{\beta k\over
2}L\right)}\left[\mbox{erf}\left({\beta kL-{\beta k_{D}\over
2}\left(\omega\tau+{1\over 2}\sin 2\omega\tau\right)\over\sqrt{2\beta
k}}\right)+\right.$ $\displaystyle\hskip 110.0pt\left.\mbox{erf}\left({\beta
kL+{\beta k_{D}\over 2}\left(\omega\tau+{1\over 2}\sin
2\omega\tau\right)\over\sqrt{2\beta k}}\right)\right]$ (3.29)

$\displaystyle I_{p}^{\rm(res)}={1\over 2~{}\mbox{erf}\left(\sqrt{\beta\over
2m}M\right)}\left[\mbox{erf}\left({{\beta\over m}M-{\beta k_{D}/m\over
4\omega}\left(1-\cos 2\omega\tau\right)\over\sqrt{2\beta/m}}\right)+\right.$
$\displaystyle\hskip 110.0pt\left.\mbox{erf}\left({{\beta\over m}M+{\beta
k_{D}/m\over 4\omega}\left(1-\cos
2\omega\tau\right)\over\sqrt{2\beta/m}}\right)\right]$ (3.30)

Because $\omega$ can be considered an intrinsic property of the system coupled
to the driving mechanism, we take it as fixed. Then, Equation (3.27), together
with (3.28)-(3.30), shows that the average of the exponential of the Jarzynski
work for a bounded ensemble of initial states depends on the protocol time
$\tau$. Indeed, for a sinusoidal protocol, there is an infinite set of values
of $\tau$ that yields the same final value $\lambda(\tau)=B$. In particular,
Eq.(3.29), shows that $I_{x}^{\rm(res)}$ may even approach $0$ or $2$, however
large $L$ is taken, for sufficiently large $\tau$. Indeed, the first error
function in Eq.(3.29) tends to $-1$, while the other tends to 1, as $\tau$
grows, all the other parameters being fixed. On the other hand, small $\tau$
implies a sum of two equal quantities, which approaches 2 for large $L$.
Tuning the values of $\tau$ one observes quite a sensitive protocol
dependence, for the average (3.27). This is illustrated in Figs. 3.2 and 3.3.
The cause of this behaviour, in presence of a resonance, is the fact that the
amplitude of the oscillator position grows linearly in time, yielding the
$\omega\tau$ term in the arguments of the error functions of
$I_{x}^{\rm(res)}$.

###  3.3 Coupled Oscillators with Periodic Forcing

In this subsection, a single oscillator, $S$, is harmonically tied to the
origin of the real line, and is harmonically driven, as in Sec. 3.2. In
addition, S is harmonically coupled to a second oscillator, E. We denote by
$k_{I}$ the stiffness of the harmonic potential linking S and E, and we also
assume that E is harmonically bound to the origin of the line, with elastic
constant $k_{E}$. Let the oscillators masses be $m_{E}$ and $m_{S}$, and let
the phase of S+E be denoted by
$\Gamma=\left(x_{E},x_{S},p_{E},p_{S}\right)=\left(\mathbf{x},\mathbf{v}\right)$,
where $p_{E}=m_{E}v_{E}$ and $p_{S}=m_{S}v_{S}$ are the momenta associated to
each oscillator. Then, the Hamiltonian of the total system is given by:

$\displaystyle{\cal H}(\mathbf{x},\mathbf{v};\lambda)$ $\displaystyle=$
$\displaystyle\left[{p_{S}^{2}\over 2m_{S}}+{k_{S}\over
2}x_{S}^{2}+{k_{D}\over 2}\left(\lambda-
x_{S}\right)^{2}\right]+\left[{p_{E}^{2}\over 2m_{E}}+{k_{E}\over
2}x_{E}^{2}\right]+{k_{I}\over 2}(x_{E}-x_{S})^{2}$ (3.31) $\displaystyle=$
$\displaystyle
H_{S}(x_{S},p_{S};\lambda)+H_{E}(x_{E},p_{E})+h_{I}(\mathbf{x})$ (3.32)

where the square brackets delimit the different contributions to the full
Hamiltonian, respectively $H_{S},\;H_{E}$ and $h_{I}$, as in Eq.(2.1) for the
JE theory. As driving term, we take the periodic protocol used above:
$\sin\gamma t$, and we set again $k=k_{S}+k_{D}$. The equations of motion for
this system are the following:

$\displaystyle\begin{cases}m_{S}\ddot{x}_{S}=-kx_{S}+k_{D}\sin\gamma
t-k_{I}\left(x_{S}-x_{E}\right)~{},\\\
m_{E}\ddot{x}_{E}=-k_{E}x_{E}-k_{I}\left(x_{E}-x_{S}\right)~{},\end{cases}$
(3.33)

with initial conditions
$\left(\mathbf{x}(0),\mathbf{\dot{x}}(0)\right)=\left(\mathbf{x}_{0},\mathbf{v}_{0}\right)=\Gamma_{0}$.
While analytical solutions for this set of equations are conceptually trivial,
they are practically involved if $k_{S}\neq k_{E}$ and $m_{S}\neq m_{E}$,
especially when integrated to compute the left hand side of the JE. On the
other hand, they can be quite simply handled in numerical calculations. We
have thus numerically sampled the initial conditions $\Gamma_{0}$ from the
truncated canonical distribution, and for each of them we have computed the
initial energy ${\cal H}(\Gamma_{0};\lambda(0))$. Then, we have numerically
solved Eqs. (3.33) for that $\Gamma_{0}$, obtaining the final condition
$\Gamma_{\tau}(\Gamma_{0})$, that has been introduced in the final Hamiltonian
${\cal H}(\Gamma_{\tau}(\Gamma_{0});\lambda(\tau))$, to obtain the work as:

$W_{J}(\Gamma_{0})={\cal H}(\Gamma_{\tau}(\Gamma_{0});\lambda(\tau))-{\cal
H}(\Gamma_{0};\lambda(0))$ (3.34)

as in Eq. (2.3), where $\lambda(0)=A=0$ and $\lambda(\tau)=\sin\gamma\tau=B$.
Collecting many works, with $\tau$ fixed, we have eventually estimated the
quantity

$\Big{\langle}e^{-\beta W_{J}}\Big{\rangle}_{0;L,M}$ (3.35)

###  3.4 Results

Our first observation is that finite size effects make the protocol dependent
the quantity (3.35) unlike the case of systems initially in contact with truly
infinite reservoirs. Of course, no real reservoir is infinite, but considering
it infinite introduces no errors when taking equilibrium averages of standard
observables, such as power laws. The situation changes if exponentials of
standard observables are considered. For the single oscillator driven by a
harmonic trap moving with constant velocity, Fig. 3.1 shows the dependence of
(3.35) on $\ell$ and on $\beta£$, for different values of the harmonic
potential stiffness $k_{p}$: fast and slow protocols yield different ensemble
averages. The cases with $L=1$ and $M=1$, represented by solid lines, show an
abrupt transition at about $\ell=1$, for small $k_{p}$. For large $k_{p}$,
dominating the coupling with the driving agent, the result gradually turns
independent of the speed of the process. Increasing the reservoir size to
$L=5$ and $M=5$, the quantity (3.35) does not appear to depend anymore on the
speed of the process $\lambda(t)$, cf. dashed lines in Fig. 3.1. In reality,
the dash-dotted lines for $L=M=2$ reveal that the process dependence merely
shifts with $L$ and $M$, becoming evident at larger $\ell$. Therefore, process
independence for (3.35) is only obtained when $L=\infty$ and $M=\infty$. An
analogous behaviour is observed as a function of the inverse temperature
$\beta$, with more evident transitions at higher temperatures.

Figure 3.1: Values of $\left\langle e^{-\beta W_{J}}\right\rangle_{L,M}$ for
the single harmonic oscillator driven by a constant speed moving harmonic
trap, with final protocol value of $\lambda(\tau)=B=1$. The result is shown as
a function of $\ell$, for different values of $L$ and $M$. Solid lines refer
to $L=M=1$, dash-dotted lines to $L=M=2$, and dashed lines to $L=M=5$. In the
left panel, blue, red, yellow and purple plots refer to
$k_{p}=0.1,\;1,\;10,\;100$, respectively. In the right panel, blue, red,
yellow and purple plots refer to $\beta=1,\;4,\;7,\;10$, respectively. Other
parameters are set to $m=1$ and $k_{D}=1$.

The second model analyzed above is even more intriguing, as resonances
significantly affect the work done on the system by external perturbations,
when finite size effects play a role. Figures 3.2 and 3.3 show that the
extension of the phase space volume does not suffice to tame the resonances
produced over sufficiently long times $\tau$. Unlike the case of infinitely
large baths, which yield the same result for all $\tau$, here a protocol
dependence arises. The reason is that a harmonic oscillator subject to no
friction and to periodic forces performs oscillations whose amplitude grows
linearly in time, if the forcing frequency equates the natural frequency of
the system. In our example, this happens for $\gamma=\omega$. Thus, the work
done on the system grows together with the amplitude, pushing $\left\langle
e^{-\beta W_{J}}\right\rangle_{L,M}$ toward $0$, at and near the resonance.

Figure 3.2: Values of $\left\langle e^{-\beta W_{J}}\right\rangle_{L,M}$ for
the single harmonic oscillator with $\lambda=\sin\gamma t$, as a function of
the forcing frequency $\gamma$, for different values of $\Gamma_{0}$ volumes
and different final times $\tau$. Left and right panels refer to $L=M=1$ and
$L=M=10$ respectively, with $\tau$ such that $B=\sin(2\pi)$ for the first
case, and $B=\sin(200\pi)$ for the second. Blue, red, yellow and purple plots
refer to inverse temperatures $\beta=1,\;10,\;100,\;1000$, respectively. Other
parameters are set to $m=1$, $k=1$, $k_{D}=1$.

Figure 3.3: Values of $\left\langle e^{-\beta W_{J}}\right\rangle_{L,M}$ for
the single harmonic oscillator with $\lambda=\sin\gamma t$, as a function of
the forcing frequency $\gamma$, for different values of $\Gamma_{0}$ volumes
and different final times $\tau$. Left and right panels refer to $L=M=1$ and
$L=M=10$ respectively, with $\tau$ such that $B=\sin(2\pi)$ for the first case
and $B=\sin(200\pi)$ for the second. Blue, red, yellow and purple plots refer
to stiffnesses $k_{s}=0.1,\;1,\;10,\;100$, respectively. Other parameters are
set to $m=1$, $\beta=1$ and $k_{D}=1$.

The inverse temperature $\beta$ and the global stiffness $k$ of the potential
also have noticeable effects on the behavior of (3.35). In particular,
increasing $\beta$ (reducing the temperature) flattens the curve about the
resonance, widening the interval of $\gamma$ values that make $\left\langle
e^{-\beta W_{J}}\right\rangle_{0;L,M}$ equal 0, rather than the infinite size
theoretical value 1, cf. Fig. 3.2. A larger value of $k$ seems instead to
stabilize $\left\langle e^{-\beta W_{J}}\right\rangle_{0;L,M}$ and reduce its
dependence on $\gamma$, as shown by Fig. 3.3.

Figure 3.4: Behaviour of $\left\langle e^{-\beta W_{J}}\right\rangle_{L,M}$ as
a function of the forcing frequency $\gamma$ for the coupled oscillators
model, simulated with $\tau$ such that $B=\sin(2\pi)$. Left, center and right
panel report the system behavior at varying values of the parameters $\beta$,
$k_{S}$ and $k_{I}$, respectively. Markers represent results from numerical
simulations, while dotted lines connecting them are a guide for the eye. Left
panel: dark blue, light blue, light grey and dark grey lines correspond to
$\beta=0.1,\;1,\;10,\;100$, respectively; other parameters are set to
$m_{S}=m_{E}=1$, $k_{S}=k_{E}=1$, $k_{I}=1$ and $L=M=1$. Center panel: dark
blue, light blue, light grey and dark grey lines correspond to
$k_{S}=0.1,\;1,\;10,\;100$, respectively; other parameters are set to
$m_{S}=m_{E}=1$, $k_{E}=1$, $k_{I}=1$, $\beta=1$ and $L=M=1$. Right panel,
dark blue and light blue lines refer to $k_{I}=10$ and $k_{I}=0.1$,
respectively, with $L=M=1$; dark grey and light grey plots refer to $k_{I}=10$
and $k_{I}=0.1$, respectively, but for $L=M=5$; remaining parameters are set
to $m_{E}=10$, $m_{S}=1$, $k_{S}=k_{E}=1$, $\beta=1$.

The pair of oscillators S and E from Sec. 3.3, with periodic forcing on S,
shows a similar behavior, at least when the coupled particles have the same
mass, $m_{S}=m_{E}=1$, and the interaction stiffness is sufficiently low (like
$k_{I}=1$). This is illustrated by the first two panels of Fig. 3.4, where the
only parameters varied are $\beta$ and $k_{S}$. It is interesting to note how
the temperature of the system needs to decrease (thus $\beta$ to increase) to
make $\left\langle e^{-\beta W_{J}}\right\rangle_{0;L,M}$ vanish. This is
especially evident in central panel of Fig. 3.4 where none of the tested
values of $k_{S}$ yields $0$ for $\beta=1$. On the contrary, $\beta=100$
obtains $\left\langle e^{-\beta W_{J}}\right\rangle_{L,M}=0$ for $k_{S}=1$ and
different values of $\gamma$. The reason is that higher $\beta$ implies a
narrower distribution, hence analogous to a case with smaller $L$ and $M$.

Note that this is relevant also for infinite baths. A small temperature causes
kinds of finite size effects, due to the smallness of the distribution
variance. That, assuming the infinite space can at least in principle be
explored, may be eliminated only at the cost of collecting enormous
statistics, which is often impossible. Therefore, the finite ensemble result
remains the only physically relevant.

The right panel of Fig. 3.4 shows the quantity (3.35) computed on a variation
of the coupled oscillators model, in which the “environment” mass $m_{E}$ is
ten times bigger than the “system” mass $m_{S}$. To include possible effects
due to the efficiency of the energy exchange between system and environment,
the stiffness of the interaction potential is varied: first a value of
$k_{I}=0.1$ is implemented, then $k_{I}=10$ is employed to account for a rapid
exchange of energy between the parts, that result almost rigidly connected.
The figures show that the two configurations generate similar results when the
initial ensemble is restricted to $L=M=1$, while noticeably different behavior
is observed for a larger system, where $L=M=5$. The rigidly coupled system
produces oscillations of $\left\langle e^{-\beta W_{J}}\right\rangle_{0;L,M}$
around the resonance frequency, with the loosely connected case exhibiting
even more evident down-peaks in the values of $\left\langle e^{-\beta
W_{J}}\right\rangle_{L,M}$ about the resonance frequency, as in the case of
the periodically forced single oscillator. Moreover, both $L=5=M$ lead to
peaks that exceed 1. Both weakly and strongly coupled oscillators indicate
that the presence of a massive environment drastically magnifies the finite
size effects, noticeably deviating from the equality $\left\langle e^{-\beta
W_{J}}\right\rangle_{0;L,M}=1$.

##  4 Irreversible Expansion of an Ideal Gas

Consider a set of $N$ identical point particles in a 2D rectangular box of
length $2L$. The particles move in straight lines, and collide elastically
with the hard boundaries of the box. The box is subdivided in two equal parts
by a wall perpendicular to two sides, that can can be removed and placed back,
according to prescribed protocols. We perform a cycle, starting from the gas
confined by the wall in the left half of the box, and in equilibrium at a
given temperature $T$. At time $t=0$, the system is isolated from the bath and
the wall is removed for a certain amount of time $\tau$. Finally, the wall is
placed back. A schematic representation of the system and its dynamics is
given in Fig. 4.5. The main observable here is the fraction of particles
trapped in the right half of the box at the end of the cycle. Because this
fraction depends on the details of the process through which the intermediate
wall is removed and reinserted, the final equilibrium state may differ from
the initial one, and may depend on the protocol.

Figure 4.5: Schematic representation of the irreversible expansion of an ideal
gas. A 2D box contains non-interacting particles in equilibrium with a heat
bath at temperature $T$. Panel (a): state of the system before the central
wall is removed; all the particles are confined in the left half of the box
and undergo specular reflections with the container walls. Panel (b): dynamics
of the particles once the central wall is removed. Panel (c): after the wall
is reintroduced, a fraction of particles is trapped in the right half of the
box.

For sake of simplicity, and without any loss of significance, because the
particles do not interact, we may replace the 2D container with a straight
line segment of length $2L$. We also assume the particles to start with random
initial positions, uniformly distributed in the interval $(-L,0)$, and with
initial velocities normally distributed, with mean $\mu=0$ and standard
deviation $\sigma=\sqrt{{k_{B}T}/{m}}$, where $k_{B}$ is the Boltzmann
constant, $T$ is the temperature of the bath and $m$ is the mass of the
particles. The fraction of particles escaping from the left half of the box to
the right half, in the time interval $[0,\tau]$, is obtained integrating over
all initial positions $x_{0}$ and initial velocities $v_{0}$ the probability
for a particle to move from $(-L,0)$ to $(0,L)$. In the limit of many
particles, the number of those leaving the left half and reaching the right
half is this probability multiplied by their total number $N$. We denote by
$N_{L}$ the number of particles in the left half of the box, and by $N_{R}$
the number of those in the right half of the box, so that $N=N_{L}+N_{R}$.
Now, note that the initial velocities pointing rightward (i.e. $v_{0}>0$) that
make a particle starting at $x_{0}\in(-L,0)$ end in $x_{\tau}\in(0,L)$ after a
time $\tau$, fulfill the inequalities:

${1\over\tau}(4nL-x_{0})<v_{0}<{1\over\tau}[(4n+2)L-x_{0}]\,,\quad
n=0,1,2,\dots$ (4.36)

because travelling a distance of $(4nL-x_{0})$ brings the particle in the
interval $(0,L)$, after a number $n$ of bounces against the left wall of the
container. Going beyond $[(4n+2)L-x_{0}]$ brings the same particle back to the
left half of the box. The same reasoning, applied to particles initially
pointing leftward (hence $v_{0}<0$), shows that a particle is found in the
right half of the box at time $\tau$ if its velocity $v_{0}$ is such that:

$-{1\over\tau}[(4n+4)L+x_{0}]<v_{0}<-{1\over\tau}[(4n+2)L+x_{0})\,,\quad
n=0,1,2,\dots$ (4.37)

The number of particles residing in the right half of the box at time
$\tau>0$, denoted as $N_{R}(\tau)$, is thus obtained by integrating over all
initial positions $x_{0}$ in the interval $(-L,0)$ and over all initial
velocities $v_{0}$ contained in the intervals given in (4.36) and (4.37). By
doing this, one implicitly assumes the system is made of infinitely many
particles, therefore the result applies only for large $N$. For sufficiently
large $N$, the following:

$\frac{N_{R}(\tau)}{N}={1\over{L\sigma\sqrt{2\pi}}}\int_{-L}^{0}dx_{0}\left\\{\sum_{n=0}^{\infty}\left[\int_{{1\over\tau}(4nL-
x_{0})}^{{1\over\tau}[(4n+2)L-x_{0}]}e^{-v_{0}^{2}/2\sigma^{2}}\mbox{d}v_{0}+\int_{{1\over\tau}[(4n+2)L+x_{0}]}^{{1\over\tau}[(4n+4)L+x_{0}]}e^{-v_{0}^{2}/2\sigma^{2}}\mbox{d}v_{0}\right]\right\\}$
(4.38)

is thus an accurate prediction for observations. Then, integrating first over
the velocity space, exchanging the integral over positions with the infinite
sum (which is made possible since each term of the infinite sum is a
continuous, bounded function), and finally integrating over initial positions,
Eq. (4.38) yields

$\displaystyle\begin{split}\frac{N_{R}(\tau)}{N}&=&{1\over
2L}\sum_{n=0}^{\infty}\left\\{\sqrt{\frac{2}{\pi}}\sigma\tau\left[\exp\left(-{16L^{2}(n+1)^{2}\over
2\sigma^{2}\tau^{2}}\right)+\exp\left(-{16L^{2}n^{2}\over
2\sigma^{2}\tau^{2}}\right)-2\exp\left(-{L^{2}(4n+2)^{2}\over
2\sigma^{2}\tau^{2}}\right)\right]\right.\\\\[5.0pt] &&\left.\hskip
10.0pt-2(4n+2)L\operatorname{erf}\left[{(4n+2)L\over\sqrt{2}\sigma\tau}\right]+4nL\operatorname{erf}\left[{4nL\over\sqrt{2}\sigma\tau}\right]+(4n+4)L\operatorname{erf}\left[{(4n+4)L\over\sqrt{2}\sigma\tau}\right]\right\\}\end{split}$
(4.39)

A comparison between this analytical result and numerical simulations of the
system is shown in Fig. 4.6. To estimate the relaxation times, in the absence
of the intermediate wall, one cannot count on the mean time between two
consecutive collisions with the boundaries of the box, given by $2L\langle
1/v\rangle$, because such a mean does not exist in a 1-dimensional space.
However, one may take the distance $2L$ divided by the mean speed
$\langle|v|\rangle$

$\widehat{t}={2L\over\langle|v|\rangle}=\sqrt{2\pi}{L\over\sigma}$ (4.40)

as a characteristic time for the dynamics. This quantity estimates the time
scale of the relaxation to a uniform distribution of particles, when their
number is sufficiently high that recurrence times can be considered infinite
to all effects. As indicated by the vertical lines in the left panel of Fig.
4.6, systems with larger $\beta$, i.e. smaller temperature, take longer to
reach the uniform distribution of particles in the box, as expected.

Figure 4.6: Left panel: behavior of $N_{R}/N$ as a function of the protocol
duration $\tau$, for different values of $\beta$, with $N=10000$, $L=5$,
$m=1$. Blue, orange and yellow solid lines refer to Eq.(4.39) for
$\beta=1,\;10,\;100$, respectively, and are obtained by truncating the formula
(4.39) to $n=1000$. The black dotted lines denote the results of the numerical
simulations. Dash-dotted vertical lines indicate the characteristic times
obtained from Eq. (4.40). Right panel: behavior of $N_{R}/N$ vs. $\tau$, for
different values of $N$, obtained from numerical simulations. Black, blue and
orange solid lines refer to $N=20,\;100,\;500$ respectively. The other
parameters are fixed to $L=5$, $\beta=10$, $m=1$.

The point here is that an irreversible process is generated by the motion of
the moving wall, that returns to its initial position, at the end of a cycle.
The result is that variations of the process time $\tau$ lead to different
results for $N_{R}$, hence of the free energy of the final equilibrium states.
For sufficiently large $\tau$, the process uniquely leads to $N_{L}=N_{R}$
(although fluctuations occur at any finite $N$, cf. the right panel of Fig.
4.6), but that is different from the initial state $N_{L}=N$, $N_{R}=0$.
Therefore, even accepting an ideal infinite thermostat, the variation of free
energy between initial and final state does not vanish and depends on the
process time $\tau$, at odds with the JE, i.e. with the canonical ensemble
from which the JE is derived. In this case, the canonical distribution of
momenta is given by a Gaussian. Were the range of momenta finite, further
corrections to the canonical prediction would arise.

For the dynamics to be Hamiltonian, as the derivation of the JE requires, the
wall could be modelled by a repulsive potential $\Phi$, placed in the center
of the box, that diverges at time 0 and $\tau$, and that vanishes at time
$\tau/2$. For instance, the following would do:

$\Phi(x;\lambda)=\left\\{\begin{array}[]{lr}0&\mbox{if
}~{}x\notin\left[1/2-\epsilon,{1}/{2}+\epsilon\right]\\\
\cfrac{1}{\lambda}-\cfrac{4}{\tau^{2}}&\mbox{if
}~{}x\in\left[1/2-\epsilon,{1}/{2}+\epsilon\right]\end{array}\right.\quad\mbox{with
}~{}~{}\lambda(t)=t(\tau-t)$ (4.41)

with $2\epsilon>0$ representing the thickness of the wall. In this case, the
infinitesimal contribution to $W_{J}$, given by the interaction of particle
$i$ in position $x_{i}$ with the wall, for a time ${\rm d}t$, is given by

${\rm d}w_{i}=\left\\{\begin{array}[]{cr}0&\mbox{if
}~{}x_{i}\notin\left[1/2-\epsilon,{1}/{2}+\epsilon\right]\\\\[4.0pt]
\cfrac{2t-\tau}{t^{2}(\tau-t)^{2}}~{}{\rm d}t&\mbox{if
}~{}x_{i}\in\left[1/2-\epsilon,{1}/{2}+\epsilon\right]\end{array}\right.$
(4.42)

which has to be integrated over the time intervals within $[0,\tau]$ such that
$x_{i}\in\left[1/2-\epsilon,{1}/{2}+\epsilon\right]$, and summed over all
particles. Now, the standard canonical formalism is not applicable to this
simple example, because the phase space corresponding to the initial
equilibrium state is altered when the wall potential is lowered to finite
height. It switches from representing an equilibrium state in half the volume
of the box, to a different equilibrium state, that occupies the whole box. The
instant in which the particles are allowed to move in whole the box, but are
still confined in its left half, the initial canonical distribution does not
describe anymore their state, and cannot be used to compute the free energy
difference between the equilibrium states before and after the wall is
removed. Nevertheless, this difficulty can be overcome, without affecting the
result, considering, as generally and efficiently done (see e.g. Section 1.3
of Ref. [27]), that physically relevant space and time scales are finite,
although they can be taken as large as one wants. Then, one realizes that a
finite but high barrier confines a finite number of particles initially in the
left half of the box, with only a negligible fraction $\epsilon$ of them
moving to the other half, for a given time. A finite, but higher barrier
confines the particles with same tolerance $\epsilon$ for a larger time, or
for the same time and a smaller $\epsilon$. Given the (arbitrarily large) time
and the (arbitrarily small) tolerance considered physically satisfactory,
there is a barrier height that produces better confinement, allowing the
initial state to be considered an equilibrium state. Then, the protocol
dependence of the free energy difference described above remains.

##  5 Concluding remarks

In this work we have discussed simple examples concerning finite size effects
and a broken time reversal symmetry on the calculations of values of
observables within the statistical mechanics formalism. It is indeed ever more
important to properly describe systems that do not lie within the traditional
bounds of statistical mechanics, developed for macroscopic systems at
equilibrium, or slowly evolving near equilibrium states. Indeed, contemporary
research widely focuses on small and far from equilibrium systems. In the case
of equilibrium macroscopic systems, the use of the standard ensembles is fully
justified, because the observables of interest are determined by the bulk, not
the tails of the probability distributions. This approach is validated both by
theory and experiments. On the contrary, fluctuations of properties of
interest in the case of small or strongly nonequilibrium systems, compare to
average signals, and require a proper characterization of the tails of the
distributions, which may also be affected by a spontaneous time reversal
symmetry breaking, not evident in the equations of motion. Indeed, the
interaction with heat reservoirs is often limited to processes that last very
short times, making effective only small parts of such environments. To
illustrate these points we have investigated simple driven harmonic
oscillators systems, averaging the popular quantity $\exp(-\beta W_{J})$ with
canonical and truncated canonical averages. We have shown that:

  * •

A single oscillator pulled by a constant speed harmonic trap yields the
infnite bath result if the process is not too fast. It sensibly and rapidly
departs from that when the speed of the driving agent grows. The effect is
more evident (as expected) for small than for large integration bounds, $L$
and $M$, for smaller harmonic constant, and for smaller bath temperatures. In
the infinite $L,M$ limits, the standard canonical result is recovered, but
larger and larger $L$ and $M$ are required, the smaller are $k_{p}$ or
$\beta$.

  * •

For a single periodically driven oscillator, the infinite bath result over a
multiple of the driving period equals 1. Strong deviations from this value,
that even reach $0$, are found for a finite bath in finite intervals about the
resonance frequency. While the theoretical result is again obtained in the
infinite $L,M$ limit, this is harder if the driving acts for longer times
(i.e., for a larger number of periods).

  * •

In the case in which the oscillator S is coupled to an oscillator E, the
infinite bath value 1 is obtained, apart from oscillations, for sufficiently
large driving frequency. Noticeable deviations from that results are still
present about specific values of the forcing frequency. In this case, we have
no analytical expression for the finite bath result. Therefore this conclusion
is based on numerical data, for a finite ensemble of initial conditions,
proven robust against variations of ensemble size.

The last example we have considered implies the breaking of the classical time
reversal symmetry and consists of an ideal gas initially in equilibrium with
an infinite thermal bath at temperature $T$, which is confined in the left
half of its container by a moving barrier. Initially, the barrier confines all
particles to the left half of the container, and that allows an equilibrium
state, represented by a uniform distribution for the positions of particles in
$[-L,0]$ and a Gaussian distribution for their velocities. As soon as the
barrier allows particles to reach the right half of the container the phase
space changes to one in which positions cover the $[-L,L]$ interval, and the
initial ensemble immediately fails to represent an equilibrium state. The
observables take instead some time to change. This prevents the application of
the standard statistical mechanical formalism, because the phase space of the
equilibrium initial state is not the one of the nonequilibrium evolution and,
for instance, the Liouville equation fails. A modification in time of the
volume of a given system to the very least requires a suitable time dependent
scaling of the phase space coordinates [28], for the formalism to apply, but
that is not possible in our case, because the volume changes instantaneously.
Nevertheless, the experiment can be performed, and a proper formalism for it
has been identified, in terms of finite potential barriers and time scales.

Discrepancies between the canonical formalism and experimental situations are
known to arise when irreversibility emerges: they are intrinsic and not merely
due to insufficient statistics. All boils down to the conclusion that finite
size and irreversibility effects similarly lead to protocol dependent averages
of exponential quantities such as $\exp(-\beta W_{J})$. The standard
statistical mechanics formalism should be adapted to treat these cases. This
fits nicely with the standard statistical mechanical justification of
ensembles. For instance, in Ref. [24], Fermi states:Studying the
thermodynamical state of a homogeneous fluid of given volume at given
temperature […] we observe that there is an infinite number of states of
molecular motion that correspond to it. With increasing time, the system
exists successively in all the dynamical states that correspond to the given
thermodynamical state. From this point of view we may say that a
thermodynamical state is the ensemble of all the dynamical states through
which, as a result of the molecular motion, the system is rapidly passing. and
Callen adds If the transition mechanism among the atomic states is
sufficiently effective, the system passes rapidly through all representative
atomic states in the course of a macroscopic observation […]. However, under
certain unique conditions, the mechanism of atomic transition may be
ineffective and the sytem may be trapped in a small subset of atypical atomic
states. Or, even if the system is not completely trapped the rate of
transition may be so slow that a macroscopic measurement does not yield a
proper average over all possible atomic states.  In reality, less than
required by Fermi and Callen is needed, for ensembles to work, because
observables of interest are generally a few and well behaved [26, 29, 30]. But
when the standard conditions are severely violated, and the observables call
for an accurate representation of large fluctuations, canonical results must
be taken with a grain of salt.

Acknowledgements. LR acknowledges financial support by the Ministero
dell’Universitá e della ricerca (Italy), grant Dipartimenti di Eccellenza 2018
– 2022 (E11G18000350001). This research was performed under the auspices of
Italian National Group of Mathematical Physics (GNFM) of INdAM.

##  Appendix A Some explicit calculations for Sec. 3.1

Retaining the notation of Sec. 2, we denote by $P_{0}(\Gamma)$ the canonical
distribution referring to a specific configuration $\Gamma$ (system $+$
environment),

$P_{0}(\Gamma)={1\over Z_{0}}e^{-\beta{\cal
H}(\Gamma;A)}\,,\quad\beta=\cfrac{1}{k_{{}_{B}}T}\;,$ (1.43)

with $k_{B}$ the Boltzmann’s constant. One thus readily finds:

$\displaystyle\hskip-55.0pt\Big{\langle}e^{-\beta
W_{J,\ell}}\Big{\rangle}_{0}={1\over
Z_{0}}\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}\mbox{d}x_{0}\mbox{d}p_{0}e^{-\beta
W_{J}(\ell;x_{0},v_{0})}e^{-\beta{\cal H}(x_{0},v_{0};0)}=$ (1.44)
$\displaystyle=\cfrac{1}{Z_{0}}\exp\left\\{-\beta\frac{k_{D}k_{p}B^{2}}{2k}-\beta{k_{D}^{2}\ell^{2}m\over
k^{2}}\left(1-\cos\omega B/\ell\right)\right\\}\times$ $\displaystyle\hskip
10.0pt\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}\mbox{d}x_{0}\mbox{d}p_{0}\exp\beta\left[{x_{0}k_{D}\ell\over\omega}\sin\omega{B\over\ell}-{p_{0}k_{D}\ell\over
k}\left(\cos\omega{B\over\ell}-1\right)-{p_{0}^{2}\over 2m}-{kx_{0}^{2}\over
2}\right]$ (1.45)

where $Z_{0}$ is the partition function of the initial canonical distribution
and it is given by

$Z_{0}=\int_{-\infty}^{\infty}\mbox{d}x\,e^{-\beta{kx^{2}/2}}\int_{-\infty}^{\infty}\mbox{d}pe^{-\beta{p^{2}/2m}}={2\pi\over\beta\omega}.$
(1.46)

The double integral in (1.45) can be separated and computed in two parts:

$\int_{-\infty}^{\infty}\mbox{d}\,x\,e^{\beta
x{k_{D}\ell\over\omega}\sin\omega{B\over\ell}-{\beta k\over
2}x^{2}}=\sqrt{2\pi\over\beta k}\exp\left({\beta k_{D}^{2}\ell^{2}m\over
2k^{2}}\sin^{2}\omega{B\over\ell}\right)$ (1.47)

and

$\int_{-\infty}^{\infty}\mbox{d}\,p\,e^{-\beta\left[p{k_{D}\ell\over
k}\left(\cos\omega{B\over\ell}-1\right)+{p^{2}\over 2m}\right]}=\sqrt{2\pi
m\over\beta}\exp\left[{\beta k_{D}^{2}\ell^{2}m\over
2k^{2}}\left(\cos^{2}\omega{B\over\ell}-2\cos\omega{B\over\ell}+1\right)\right]\;,$
(1.48)

from which it follows that:

$\displaystyle\hskip-50.0pt\Big{\langle}e^{-\beta
W_{J,\ell}}\Big{\rangle}_{0}={\beta\omega\over
2\pi}{2\pi\over\beta\omega}\exp\left\\{-\beta{k_{D}k_{p}B^{2}\over
2k}-\beta{k_{D}^{2}\ell^{2}\over k\omega}\left(1-\cos\omega
B/\ell\right)\right\\}\times$ $\displaystyle\hskip 15.0pt\exp\left[{\beta
k_{D}^{2}\ell^{2}\over 2k\omega}\sin^{2}\omega{B\over\ell}+{\beta
k_{D}^{2}\ell^{2}\over
2k\omega}\left(\cos^{2}\omega{B\over\ell}-2\cos\omega{B\over\ell}+1\right)\right]$
$\displaystyle=\exp\left\\{-\beta{k_{D}k_{p}B^{2}\over 2k}\right\\}\;.$ (1.49)

Let us now turn to consider canonical distributions truncated at a given
distance $L$ from the rest position of the oscillator, and at a maximum
momentum $M$. Referring to the model of a single oscillator subject to a
linear protocol, treated in Sec. 3.1, we denote:

$P_{0}(x,p)=\cfrac{1}{Z_{0}(L,M)}\left\\{\begin{array}[]{lcr}e^{-\beta{(kx^{2}+p^{2}/m)/2}}&\mbox{if}&|x|\leq
L~{}\mbox{ and }~{}|p|\leq M\\\\[7.0pt] 0&\mbox{if}&|x|>L~{}\mbox{ or
}~{}|p|>M\end{array}\right.$ (1.50)

where:

$Z_{0}(L,M)=\int_{-L}^{L}\mbox{d}x\,e^{-\beta{kx^{2}/2}}\int_{-M}^{M}\mbox{d}pe^{-\beta{p^{2}/2m}}={2\pi\over\beta\omega}\,\mbox{erf}\left(\sqrt{\beta
k\over 2}L\right)\mbox{erf}\left(\sqrt{\beta\over 2m}M\right)$ (1.51)

Consequently, the average of $\exp(-\beta W_{J})$ for a given $\ell$ now
reads:

$\displaystyle\hskip-60.0pt\Big{\langle}e^{-\beta
W_{J,\ell}}\Big{\rangle}_{0;L,M}=\cfrac{1}{Z_{0;L,M}}\exp\left[-\beta{k_{D}k_{p}B^{2}\over
2k}-\beta{k_{D}^{2}\ell^{2}m\over k^{2}}\left(1-\cos\omega
B/\ell\right)\right]\times$ (1.52)
$\displaystyle\hskip-30.0pt\int_{-L}^{L}\int_{-M}^{M}\mbox{d}x_{0}\mbox{d}p_{0}\exp\left\\{\beta\left[{x_{0}k_{D}\ell\over\omega}\sin\omega{B\over\ell}-{p_{0}k_{D}\ell\over
k}\left(\cos\omega{B\over\ell}-1\right)-{p_{0}^{2}\over 2m}-{kx_{0}^{2}\over
2}\right]\right\\}$ (1.53)

where we can separately compute:

$\displaystyle\hskip-50.0pt\int_{-L}^{L}\mbox{d}\,x\,e^{\beta
x{k_{D}\ell\over\omega}\sin\omega{B\over\ell}-{\beta k\over
2}x^{2}}=\sqrt{\pi\over 2\beta k}e^{\beta
k_{D}^{2}\ell^{2}m\sin^{2}\omega{B\over\ell}\over 2k^{2}}$ (1.54)
$\displaystyle\times\left[\mbox{erf}\left(\sqrt{\beta k\over 2}L+\sqrt{\beta
m\over 2}{k_{D}\ell\sin\omega{B\over\ell}\over
k}\right)+\mbox{erf}\left(\sqrt{\beta k\over 2}L-\sqrt{\beta m\over
2}{k_{D}\ell\sin\omega{B\over\ell}\over k}\right)\right]$ (1.55)

and

$\displaystyle\hskip-65.0pt\int_{-M}^{M}\mbox{d}\,p\,e^{-\beta\left[p{k_{D}\ell\over
k}\left(\cos\omega{B\over\ell}-1\right)+{p^{2}\over 2m}\right]}=\sqrt{\pi
m\over 2\beta}e^{{\beta k_{D}^{2}\ell^{2}m\over
2k^{2}}\left(\cos\omega{B\over\ell}-1\right)^{2}}$ (1.56)
$\displaystyle\times\left[\mbox{erf}\left(\sqrt{\beta\over 2m}M-\sqrt{\beta
m\over 2}{k_{D}\ell\over
k}\left(\cos\omega{B\over\ell}-1\right)\right)\right.$ (1.57)
$\displaystyle\hskip 85.0pt\left.+\,\mbox{erf}\left(\sqrt{\beta\over
2m}M+\sqrt{\beta m\over 2}{k_{D}\ell\over
k}\left(\cos\omega{B\over\ell}-1\right)\right)\right]$ (1.58)

Therefore, one obtains:

$\Big{\langle}e^{-\beta W_{J,\ell}}\Big{\rangle}_{0;L,M}=I_{exp}\cdot
I_{x}\cdot I_{p}$ (1.59)

with

$\displaystyle I_{exp}=\exp\left\\{-\beta{k_{D}k_{p}B^{2}\over
2k}\right\\}=\Big{\langle}e^{-\beta W_{J,\ell}}\Big{\rangle}_{0}$ (1.60)
$\displaystyle I_{x}={\mbox{erf}\left(\sqrt{\beta k\over 2}L+\sqrt{\beta
m\over 2}{k_{D}\over
k}\ell\sin\omega{B\over\ell}\right)+\mbox{erf}\left(\sqrt{\beta k\over
2}L-\sqrt{\beta m\over 2}{k_{D}\over k}\ell\sin\omega{B\over\ell}\right)\over
2~{}\mbox{erf}\left(\sqrt{\beta k\over 2}L\right)}$ (1.61) $\displaystyle
I_{p}={\mbox{erf}\left(\sqrt{\beta\over 2m}M+\sqrt{\beta m\over 2}{k_{D}\over
k}\ell\left(\cos\omega{B\over\ell}-1\right)\right)+\,\mbox{erf}\left(\sqrt{\beta\over
2m}M-\sqrt{\beta m\over 2}{k_{D}\over
k}\ell\left(\cos\omega{B\over\ell}-1\right)\right)\over
2~{}\mbox{erf}\left(\sqrt{\beta\over 2m}M\right)}$ (1.62)

##  Appendix B Some explicit calculations for Sec. 3.2

For the model of a single oscillator subject to a periodic forcing, discussed
in Sec. 3.2, one has:

$\displaystyle\hskip-40.0pt\Big{\langle}e^{-\beta
W_{J}}\Big{\rangle}_{0;L,M}={1\over Z_{0}(L,M)}\times$
$\displaystyle\hskip-25.0pt\exp\left\\{-{\beta k_{D}\over
4}\left(1-{k_{D}/m\over\omega^{2}-\gamma^{2}}\right)\left(1-\cos
2\gamma\tau\right)-{\beta
k_{D}^{2}\gamma^{2}/m\over(\omega^{2}-\gamma^{2})^{2}}\left(1-{\gamma\over\omega}\sin\gamma\tau\sin\omega\tau-\cos\gamma\tau\cos\omega\tau\right)\right\\}\times$
$\displaystyle\int_{-L}^{L}\mbox{d}x\,\exp\left\\{-{\beta k\over
2}x^{2}+{\beta
k_{D}\gamma\omega\over\omega^{2}-\gamma^{2}}\left(\cos\gamma\tau\sin\omega\tau-{\gamma\over\omega}\sin\gamma\tau\cos\omega\tau\right)x\right\\}\times$
$\displaystyle\int_{-M}^{M}\mbox{d}p\,\exp\left\\{-{\beta\over 2m}p^{2}+{\beta
k_{D}\gamma/m\over\omega^{2}-\gamma^{2}}\left(1-{\gamma\over\omega}\sin\gamma\tau\sin\omega\tau-\cos\gamma\tau\cos\omega\tau\right)p\right\\}$
(2.63)

which thus leads to:

$\Big{\langle}e^{-\beta W_{J}}\Big{\rangle}_{0;L,M}=I_{exp}\cdot I_{x}\cdot
I_{p}$ (2.64)

where

$\displaystyle\hskip-50.0ptI_{exp}=\exp\left\\{-{\beta k_{D}\over
4}\left(1-{k_{D}/m\over\omega^{2}-\gamma^{2}}\right)\left(1-\cos
2\gamma\tau\right)\right.$ (2.65) $\displaystyle\hskip 20.0pt\left.-{\beta
k_{D}^{2}\gamma^{2}\over(\omega^{2}-\gamma^{2})^{2}}\left[{1\over
m}\left(1-{\gamma\over\omega}\sin\gamma\tau\sin\omega\tau-\cos\gamma\tau\cos\omega\tau\right)\right.\right.$
(2.66) (2.67) $\displaystyle\hskip 20.0pt\left.-{1\over
2m}\left(1-{\gamma\over\omega}\sin\gamma\tau\sin\omega\tau-\cos\gamma\tau\cos\omega\tau\right)^{2}\right.$
(2.68) $\displaystyle\hskip 40.0pt\left.-{\omega^{2}\over
2k}\left(\cos\gamma\tau\sin\omega\tau-{\gamma\over\omega}\sin\gamma\tau\cos\omega\tau\right)^{2}\biggr{]}\vphantom{\frac{1}{2}}\right\\}$
(2.69)

$\displaystyle I_{x}={1\over 2~{}\mbox{erf}\left(\sqrt{\beta k\over
2}L\right)}\left[\mbox{erf}\left({\beta kL-{\beta
k_{D}\gamma\omega\over\omega^{2}-\gamma^{2}}\left(\cos\gamma\tau\sin\omega\tau-{\gamma\over\omega}\sin\gamma\tau\cos\omega\tau\right)\over\sqrt{2\beta
k}}\right)+\right.$ (2.70) $\displaystyle\hskip
105.0pt\left.\mbox{erf}\left({\beta kL+{\beta
k_{D}\gamma\omega\over\omega^{2}-\gamma^{2}}\left(\cos\gamma\tau\sin\omega\tau-{\gamma\over\omega}\sin\gamma\tau\cos\omega\tau\right)\over\sqrt{2\beta
k}}\right)\right]$ (2.71)

$\displaystyle I_{p}={1\over 2~{}\mbox{erf}\left(\sqrt{\beta\over
2m}M\right)}\left[\mbox{erf}\left({{\beta\over m}M-{\beta
k_{D}\gamma/m\over\omega^{2}-\gamma^{2}}\left(1-{\gamma\over\omega}\sin\gamma\tau\sin\omega\tau-\cos\gamma\tau\cos\omega\tau\right)\over\sqrt{2\beta/m}}\right)\right.$
(2.72) $\displaystyle\hskip 105.0pt\left.+\mbox{erf}\left({{\beta\over
m}M+{\beta
k_{D}\gamma/m\over\omega^{2}-\gamma^{2}}\left(1-{\gamma\over\omega}\sin\gamma\tau\sin\omega\tau-\cos\gamma\tau\cos\omega\tau\right)\over\sqrt{2\beta/m}}\right)\right]\;.$
(2.73)

## References

  * [1] C. H. Bennett, “Efficient estimation of free energy differences from monte carlo data,” Journal of Computational Physics, vol. 22, no. 2, pp. 245–268, 1976.
  * [2] B. Widom, “Some topics in the theory of fluids,” The Journal of Chemical Physics, vol. 39, no. 11, pp. 2808–2812, 1963.
  * [3] R. W. Zwanzig, “High-temperature equation of state by a perturbation method. ii. polar gases,” The Journal of Chemical Physics, vol. 23, no. 10, pp. 1915–1922, 1955.
  * [4] C. Jarzynski, “Nonequilibrium equality for free energy differences,” Physical Review Letters, vol. 78, no. 14, p. 2690, 1997.
  * [5] G. E. Crooks, “Entropy production fluctuation theorem and the nonequilibrium work relation for free energy differences,” Physical Review E, vol. 60, no. 3, p. 2721, 1999.
  * [6] C. Giberti, L. Rondoni, and C. Vernia, “Temperature and correlations in 1-dimensional systems,” The European Physical Journal Special Topics, vol. 228, pp. 129–142, 2019.
  * [7] C. Giberti, L. Rondoni, and C. Vernia, “O (n) fluctuations and lattice distortions in 1-dimensional systems,” Frontiers in Physics, vol. 7, p. 180, 2019.
  * [8] G. Dematteis, L. Rondoni, D. Proment, F. De Vita, and M. Onorato, “Coexistence of ballistic and fourier regimes in the $\beta$ fermi-pasta-ulam-tsingou lattice,” Physical Review Letters, vol. 125, no. 2, p. 024101, 2020.
  * [9] M. Lax, Symmetry Principles in Solid State and Molecular Physics.  John Wiley & Sons, 1972.
  * [10] R. G. Sachs, The physics of time reversal.  University of Chicago Press, 1987.
  * [11] M. Robnik and M. V. Michael V Berry, “False time-reversal violation and energy level statistics: the role of anti-unitary symmetry,” J. Phys. A: Math. Theor., vol. 19, p. 669, 1986.
  * [12] J. Roberts and G. Quispel, “Chaos and time-reversal symmetry. order and chaos in reversible dynamical systems,” Phys. Rep., vol. 216, pp. 63–177, 1992\.
  * [13] E. G. D. Cohen and G. Gallavotti, “Note on two theorems in nonequilibrium statistical mechanics,” Journal of Statistical Physics, vol. 96, pp. 1343–1349, 1999.
  * [14] D. Evans, E. Coehn, and G. P. Morriss, “Probability of second law violations in shearing steady states,” Phys. Rev. Lett., vol. 71, p. 2401, 1993.
  * [15] M. Barbier and P. Gaspard, “Microreversibility, nonequilibrium current fluctuations, and response theory,” Journal of Physics A: Mathematical and Theoretical, vol. 51, p. 355001, 2018.
  * [16] L. Onsager, “Reciprocal relations in irreversible processes. i.,” Phys. Rev., vol. 37, no. 4, p. 405, 1931.
  * [17] L. Onsager, “Reciprocal relations in irreversible processes. ii.,” Phys. Rev., vol. 38, no. 12, p. 2265, 1931.
  * [18] H. B. Callen and T. A. Welton, “Irreversibility and generalized noise,” Phys. Rev., vol. 83, p. 34, 1951.
  * [19] R. Kubo, “Statistical-mechanical theory of irreversible processes. i. general theory and simple applications to magnetic and conduction problems,” Journal of the Physical Society of Japan, vol. 12, pp. 570–586, 1957.
  * [20] R. Kubo, “The fluctuation-dissipation theorem,” Reports on progress in physics, vol. 29, p. 255, 1966.
  * [21] H. B. G. Casimir, “On onsager’s principle of microscopic reversibility,” Reviews of Modern Physics, vol. 17, p. 343, 1945.
  * [22] J. Varillas, G. Ciccotti, J. Alcalá, and L. Rondoni, “Jarzynski equality on work and free energy: Crystal indentation as a case study,” The Journal of Chemical Physics, vol. 156, no. 11, p. 114118, 2022.
  * [23] U. Seifert, “Stochastic thermodynamics, fluctuation theorems and molecular machines,” Rep. Prog. Phys., vol. 75, no. 12, p. 126001, 2012.
  * [24] E. Fermi, Thermodynamics.  Lect. Notes Phys., Dover Books on Physics, 2012.
  * [25] H. B. Callen, Thermodynamics and an introduction to thermostatistics.  Wiley, 1985.
  * [26] A. Khinchin, Mathematical Foundations of Statistical Mechanics.  Martino Fine Books, 2014.
  * [27] R. Kubo, M. Toda, and N. Hashitsume, Statistical physics II: nonequilibrium statistical mechanics, vol. 31.  Springer Science & Business Media, 2012.
  * [28] G. Ciccotti and L. Rondoni, “Jarzynski on work and free energy relations: The case of variable volume,” AIChE Journal, vol. 67, no. 1, p. e17082, 2021\.
  * [29] U. M. B. Marconi, A. Puglisi, L. Rondoni, and A. Vulpiani, “Fluctuation–dissipation: response theory in statistical physics,” Physics reports, vol. 461, no. 4-6, pp. 111–195, 2008.
  * [30] L. Rondoni, Frontiers and Progress of Current Soft Matter Research, ch. Introduction to Nonequilibrium Statistical Physics and Its Foundations.  New York: Springer, 2021.


# Can Identifier Splitting Improve Open-Vocabulary Language Model of Code?

Jieke Shi, Zhou Yang, Junda He, Bowen Xu$\ast$, David Lo $\ast$ Corresponding
author. School of Computing and Information Systems, Singapore Management
University
{jiekeshi, zyang, jundahe, bowenxu.2017<EMAIL_ADDRESS>
###### Abstract

Statistical language models on source code have successfully assisted software
engineering tasks. However, developers can create or pick arbitrary
identifiers when writing source code. Freely chosen identifiers lead to the
notorious out-of-vocabulary (OOV) problem that negatively affects model
performance. Recently, Karampatsis et al. showed that using the Byte Pair
Encoding (BPE) algorithm to address the OOV problem can improve the language
models’ predictive performance on source code. However, a drawback of BPE is
that it cannot split the identifiers in a way that preserves the meaningful
semantics. Prior researchers also show that splitting compound identifiers
into sub-words that reflect the semantics can benefit software development
tools. These two facts motivate us to explore whether identifier splitting
techniques can be utilized to augment the BPE algorithm and boost the
performance of open-vocabulary language models considered in Karampatsis et
al.’s work.

This paper proposes to split identifiers in both constructing vocabulary and
processing model inputs procedures, thus exploiting three different settings
of applying identifier splitting to language models for the code completion
task. We contrast models’ performance under these settings and find that
simply inserting identifier splitting into the pipeline hurts the model
performance, while a hybrid strategy combining identifier splitting and the
BPE algorithm can outperform the original open-vocabulary models on predicting
identifiers by $3.68\%$ of recall and $6.32\%$ of Mean Reciprocal Rank. The
results also show that the hybrid strategy can improve the entropy of language
models by $2.02\%$.

###### Index Terms:

Open Vocabulary, Identifier Splitting, Language Model of Code

##  I Introduction

Numerous works have applied statistical language models (LMs) on source code
to help tackle important tasks in software engineering, including code
completion [1], program repair [2], and many others [3]. Same as modeling
natural language, creating appropriate vocabulary is a crucial prerequisite
[4]. However, when writing source code, software developers can create
arbitrary identifiers they like, which probably contain multiple words, e.g.
addItemsToList. Due to this fundamental fact, models of code often get an
extremely sparse vocabulary containing many rare words when processing code
corpora. Training models with such sparse (and typically large) vocabulary is
ineffective, and obtained models often have poor performance [4]. In addition,
if identifiers are not observed in the vocabulary, the model cannot handle
them, which is known as the notorious out-of-vocabulary (OOV) problem.

Currently, open-vocabulary methods like Byte-Pair Encoding (BPE) algorithm [5]
are widely used in modeling natural languages and achieve promising results in
practice. These methods can solve the OOV problem while customizing the size
of the vocabulary. Inspired by such success, Karampatsis et al. [6] first
applied the BPE algorithm to construct vocabulary from source code and showed
that open-vocabulary LMs have outstanding performance on the code completion
task. However, BPE selects the most frequent sub-words into the vocabulary,
and this frequency-based approach often fails to capture the semantics and
intentions of identifier names when choosing sub-words. Although developers
create any identifiers at will, they usually follow certain naming conventions
that make identifiers meaningful, legible and easy to understand, either in
camelCase or in snake_case [7]. For example, the method name getListener
follows the camelCase convention, and a programmer can easily infer that this
method can be used to get a Listener object. At the same time, the BPE
algorithm will represent it as three sub-words in our preliminary study: get,
List and ener, the latter two do not reflect the semantics developers try to
convey.

In order to empower the BPE algorithm with the ability to better sense
semantics when splitting words, an intuitive preprocessing strategy is to
split compound identifiers into several words that can imply certain meanings,
which is called identifier splitting techniques. Prior research works have
demonstrated that various information retrieval models for program
comprehension tasks are benefit from identifier splitting, e.g., feature-
related code localization [8], code reuse [9]. However, this empirical
conclusion is ambiguous for modern LMs of code as no such work has
demonstrated it. As stated above, open-vocabulary LMs assist many software
engineering tasks effectively but are weak in capturing the semantics of
identifiers when creating vocabulary. Thus, it is imperative to clarify if we
can improve the performance at a more considerable margin by combining
identifier splitting with the BPE algorithm.

In this paper, we investigate the potential benefits of splitting identifiers
in open-vocabulary LMs of code. Specifically, we adopt the same LMs presented
by Karampatsis et al. [6], which are the first to adopt the BPE algorithm in
code modeling. To achieve the goal, we propose to apply identifier splitting
in two stages of open-vocabulary LMs: vocabulary construction and model input
processing. Furthermore, we propose two different preprocessing strategies in
these stages to apply identifier splitting techniques: (1) simple strategy: we
split all identifiers in vocabulary construction and apply identifier
splitting before in model input processing stages; (2) hybrid strategy: in the
vocabulary construction stage, we first split all identifiers and then merge
them with original corpora. In the model input processing stage, we apply
identifier splitting only when BPE fails to tokenize them as the original
forms. We train LMs under these different settings and evaluate them on the
code completion task as [6] to show the effectiveness of identifier splitting
in open-vocabulary LMs.

We perform experiments on the C language dataset released by Karampatsis et
al. [6]. We evaluate the cross entropy of LMs, and use Mean Reciprocal Rank
(MRR) to measure the performance of LMs on the code completion task.
Furthermore, we also obtain the MRR and recall at rank 10 (R@10) on predicting
identifier tokens (excluding keywords, punctuations, etc.). The experimental
results show that simply performing identifier splitting into preprocessing
procedures does not suffice; it degrades MRR by $0.46\%$ and $5.68\%$ on
predicting all tokens and identifiers, respectively. At the same time, the
hybrid strategy is more effective for open-vocabulary LMs, outperforming the
LMs with the original setting by $6.23\%$ of MRR on predicting identifier
tokens. The improvements of $2.02\%$ of LMs entropy and $3.68\%$ in terms of
R@10 on predicting identifiers also confirm the hybrid strategy’s
effectiveness. The results highlight that the identifier splitting can be
combined with open-vocabulary methods to enhance the performance of language
models of source code.

The rest of this paper is organized as follows. Section II briefly describes
backgrounds of this paper. In Section III, we elaborate our methodology to
apply identifier splitting in the code completion pipelines. We describe the
experiment settings and present the results of our experiments in Section IV.
Section V discusses some related works. Finally, we conclude the paper and
present future work in Section VI.

##  II Background

This section introduces the background of related techniques in this work,
including identifier splitting techniques and the Byte-Pair Encoding (BPE)
algorithm.

###  II-A Naming Convention and Identifier Splitting

While identifier names do not affect program functions and human developers
can pick or create names at will, some naming conventions are encouraged to
follow as they can give meaningful and well-readable names so that other
programmers’ comprehension of code can be considerably improved [10]. Two
widely adopted naming conventions are camelCase rule and snake_case rule [7].
However, simply splitting by conventions is not accurate enough when
extracting meaningful sub-words from compound identifiers, especially for
identifiers that do not strictly follow naming conventions (e.g., the same-
case identifier where all characters are in the single case like httprequest).
Several more precise and innovative identifier splitting approaches are
proposed to handle more sophisticated situations beyond conventions, in which
Ronin [11] is the state-of-the-art. Ronin splits identifiers into sub-words
based on various heuristic rules and a pre-defined frequent sub-word table. In
this paper, we take Ronin as a representative of identifier splitting tools.

###  II-B Byte-Pair Encoding

In this paper, we use a popular open-vocabulary method called the Byte-Pair
Encoding (BPE) algorithm [5]. BPE algorithm consists of two components: (1)
the vocabulary construction stage, which takes text corpora and returns a
vocabulary with the predefined size; and (2) the tokenization stage, which
segments and tokenizes new corpora with the built vocabulary and returns a
sequence of tokens. Initially, corpora are split into a set of tokens while
each token only contains one character. The BPE algorithm iteratively merges
the most frequent pair of tokens into a new single token until a given maximum
number of merge operations is reached. These single tokens are then used to
replace the original pair of tokens in the corpora. The resulting vocabulary
is an ordered list of sub-word units created from the merge operations. When
we run the BPE algorithm on new corpora using the vocabulary, the tokens are
merged in the same order as occurred during vocabulary construction.

##  III Methodology

Figure 1: The overview of three strategies about applying the BPE algorithm
and identifier splitting in the language models. Original refers to the
original setting, which only adopts BPE, and no identifier splitting is
involved. Simple Strategy and Hybrid Strategy refer to two different
strategies which combine identifier splitting with BPE.

Figure 1 presents an overview of how we apply identifier splitting in the
open-vocabulary language models. It shows that identifier splitting can be
used in both vocabulary construction and model inputs processing procedures
before training the language models. We elaborate the motivation and
processing details as follows.

Vocabulary Construction. As shown in Figure 1, we apply the BPE algorithm on
the input code corpora and output a vocabulary with a predefined size. The
vocabulary will be used in the tokenization of model inputs later.

As introduced in Section II, the BPE algorithm builds a vocabulary based on
the frequency of pairs appearing in training corpora. Prior works sample some
software projects as a corpus to create a vocabulary, and then use the
vocabulary to tokenize other mutually exclusive projects [6]. We notice the
fact that many complex identifiers are usually project-specific or even file-
specific. Researchers have shown that code has a high degree of localness,
where identifiers are repeated often within close distance while they are
rarely used in other projects [12]. It means that the vocabulary created on
one corpus may contain many complex compounds (as they appear very frequently
in the corpus). However, those complex compounds are likely to be rare in
other corpora. For a vocabulary with limited size, these rare compounds
infringe the space of sub-words that could be applied in the BPE procedure of
other corpora, reducing the efficiency of the vocabulary. It motivates us to
apply identifier splitting techniques to the corpus.

Identifier splitting can decompose complex identifiers into several sub-words.
Karampatsis et al. [6] find that on a corpus of around $11.6$ million unique
tokens, the size of this corpus decrease dramatically by around $90\%$ after
splitting identifiers. It implies that although the complex compounds are
usually project-specific, the subunits that make up these compounds are highly
repetitive across different projects. Instead of creating a vocabulary that
contains many complex words, we propose to construct a vocabulary based on the
corpus after identifier splitting (Split Corpora as in Figure 1).

We do not claim that using identifier splitting to process the corpus is
always beneficial. Integrating identifier splitting in the vocabulary creation
may lead to some negative impacts. The vocabulary created by the corpora after
identifier splitting may need to use multiple sub-words to represent common
compounds shared across different projects, e.g., popular API names. Splitting
such common compounds may increase the length of the tokenized sequence and
make it harder to relate the current prediction to the past context of inputs
[6]. Thus, we propose another strategy that merges the original and split
corpora and builds a vocabulary from the combined corpora to tackle the issue
(Merged Corpora as in Figure 1).

Model Inputs Processing. As shown in Figure 1, the model input processing part
takes input tokens and processes them into lists of sub-tokens for language
model training. To handle an input token that is not in the vocabulary, we use
the BPE algorithm to decompose this input into a list of sub-tokens and then
feed these sub-tokens into the model. As stated in Section II, the vocabulary
created by the BPE algorithm is an ordered list of sub-words. When BPE
decomposes input tokens, it will follow the same order as recording in the
vocabulary. We take the identifier getCategory as an example. If we directly
apply BPE (using a vocabulary of $10,000$ words) to this identifier, we get
the following three sub-tokens: getC, ateg and ory, which obviously break the
original semantics of the identifier name. Although the word Category is in
vocabulary, its position ($9096th$) is almost close to the end. When BPE
traverses the vocabulary, it will encounter and create the sub-token getC
($1383th$) much earlier. If we first split this identifier into get and
Category, and then apply BPE (using the same vocabulary) on the two words, we
still get get and Category. Identifier splitting can utilize the semantic
information conveyed with naming convention and prevent less meaningful sub-
tokens (e.g., getC and ateg) from being created. This observation inspires us
to apply identifier splitting before BPE.

Splitting model inputs may also lead to negative impacts. For instance,
identifiers (e.g., types of exceptions or methods like toString) can be shared
across different projects, especially in object-oriented programming languages
like Java. Such identifiers can frequently appear in the corpus and
consequently are included in the vocabulary. They can be compactly represented
only using one token, while identifier splitting will force them to be
represented using multiple sub-tokens. As a result, we use a hybrid strategy
to mitigate such negative impacts. More specifically, we first apply BPE to an
identifier. If the tokenized result is identical to the original identifier,
we directly feed it into the model. Otherwise, we feed the separated tokens to
BPE after applying identifier splitting.

Considering the above, we combine the different operations in the vocabulary
construction and the model inputs and propose the following three settings as
shown in Figure 1 to explore the effectiveness of identifier splitting in the
open-vocabulary LMs:

  * •

Original: using the BPE algorithm to create a vocabulary directly and then use
the vocabulary to tokenize corpora as input. No identifier splitting is
applied in this setting;

  * •

Simple strategy: splitting all identifiers in corpora first then use BPE to
construct a vocabulary, and splitting all identifiers in model inputs;

  * •

Hybrid strategy: splitting identifiers and merging them with original corpora
for BPE vocabulary construction, and splitting identifiers in model inputs
only when BPE fails to tokenize them as the original forms.

##  IV Experiments and Result Analysis

###  IV-A Implementation and Datasets

To make the experiments under a computationally feasible scale, we select a
lightweight, yet still effective Gated Recurrent Unit (GRU) model [13]. The
model is also used in a recent work by Karampatsis et al. [6], which aims to
analyze how the BPE algorithm can improve LMs of code. Also, inspired by [6],
we limit the vocabulary size to 10k, set the input length and dimension of the
GRU model as $200$ and $512$, and train all models using the stochastic
gradient descent optimizer with a learning rate of $0.1$ and a mini-batch size
of $32$.

We use the dataset released by Karampatsis et al. [6] in our experiments,
which consists of $177/141/73$ open-source projects in C language for
training/validation/testing. Before training LMs, the corpora are processed in
the same way as [6], removing strings of more than $15$ characters length,
non-ASCII tokens and comments. We only use the training set to construct the
vocabularies. The replication package are available via
https://github.com/soarsmu/CodeNLM.git.

###  IV-B Target Task and Evaluation Metrics

Language Model. We use the average per token cross entropy to evaluate the
performance of our language models. The cross entropy is viewed as an
intrinsic metric of LMs and employed in the previous work [6]. By computing
the logarithm mean of probability scores assigned by the LM over a sequence of
source code, it estimates the average of bits required when using LMs to
predict each token. A lower value of the cross entropy is favourable because
it indicates LMs are easier to make correct predictions. Because the open-
vocabulary LMs are based on sub-words units, the cross entropy is formalized
as follows to compute the distribution over all sub-words
$w_{i}^{1},\dots,w_{i}^{m-1}$ instead of each token $t_{i}$:

$\small
H(N)=-\frac{1}{N}\sum_{n=1}^{N}log\prod_{m=1}^{M}p(w_{i}^{m}|t_{1},...,t_{n-1},w_{i}^{1},...,w_{i}^{m-1})$
(1)

where $N$ is the number of tokens in the sequence, $M$ is the number of sub-
words contained in $t_{i}$, and
$p(w_{i}^{m}|t_{1},...,t_{n-1},w_{i}^{1},...,w_{i}^{m-1})$ is the probability
of the sub-word $w_{i}^{m}$ given all the previous tokens and sub-words.

Code Completion. Automated code completion, an essential feature of modern
integrated development environments (IDEs), aims to suggest a range of
possible subsequent tokens within a toggle list. In open-vocabulary models,
while a complete token could be a combination of multiple sub-words, this task
can be formalized as maximizing the probability: $\arg\max
p(w_{1},\cdots,w_{n}|code\\_before)$, in which $code\\_before$ is the previous
code snippet and $w_{1},\cdots,w_{n}$ constitute the next complete token. The
probability of the next complete token is the product of the probability of
each sub-word. To obtain a complete token in open-vocabulary LMs, we use a
customized beam search algorithm introduced by [6], which can efficiently
search through the sub-word expansion space and returns the top $k$ most
possible complete candidates.

We evaluate the results with the widely-used Mean Reciprocal Rank (MRR)
metric. MRR takes the rank of the correct answer as the primary grading
criteria. For each token, if the correct answer ranks $n$th position among the
top $k$ candidates, the score would be $\frac{1}{n}$. MRR is calculated by the
following equation.

$\small MRR=\frac{1}{|T|}\sum_{i=1}^{|k|}\frac{1}{rank_{i}}$ (2)

Statistics conducted by Hellendoorn et al. [14] on real-world code completion
scenarios observes that LMs for completion perform worse on identifiers than
other types of tokens. Therefore, we also present the MRR and recall at rank
10 (R@10) on predicting identifier tokens particularly (excluding keywords,
punctuations, etc.).

TABLE I: Performance of the considered models under different strategies. Strategy | All Tokens | Identifiers
---|---|---
Entropy | MRR | R@10 | MRR
Original | 4.46 | 64.61 | 37.55 | 21.83
Simple | 4.45(-0.22%) | 64.31(-0.46%) | 36.26(-3.44%) | 20.59(-5.68%)
Hybrid | 4.37(-2.02%) | 65.24(+0.98%) | 38.93(+3.68%) | 23.19(+6.23%)

###  IV-C Results

We compare the performance of language models with different strategies. We
denote the LMs with original settings, LMs with simple strategy and LMs with
hybrid strategy as OriLMs, SSLMs and HSLMs, respectively.

Table I shows the performance of different LMs. We find that SSLMs perform the
worst among the three models. Although the entropy of SSLMs is slightly
improved in comparison to OriLMs, SSLMs degrade the MRR on predicting all
tokens by $0.46\%$, which reflects that the model’s performance is not good as
the OriLMs. For predicting identifier tokens, we observe a bigger performance
gap between OriLMs and SSLMs. In terms of R@10 and MRR on predicting
identifiers, OriLMs outperform SSLMs by $3.44\%$ and $5.68\%$, respectively.
These results indicate that the simple strategy is not adequate for open-
vocabulary LMs and even hurts the entropy and performance of LMs on the code
completion task in most cases.

However, we observe that HSLMs outperform OriLMs by up to $0.98\%$ in terms of
MRR on predicting all tokens. The models’ entropy decreases $2.02\%$, showing
that the hybrid strategy boosts open-vocabulary LMs. Compared with the results
on all token prediction, HSLMs outperform OriLMs by a larger margin on
identifier prediction: the results are boosted to $3.68\%$ and $6.23\%$ in
terms of R@10 and MRR, respectively. The results demonstrate that combining
identifier splitting with the BPE algorithm can improve the performance of
open-vocabulary LMs, especially the improved performance on identifier
prediction reveals that the LMs with hybrid strategy can synthesize
identifiers from sub-words better.

In summary, by following the hybrid strategy, identifier splitting can boost
the performance of open-vocabulary LMs of code. The performance of open-
vocabulary LMs can be improved by $0.98\%$ and $6.23\%$ in terms of MRR on
predicting both all tokens and identifiers. The entropy of LMs and R@10 on
predicting identifiers can also be improved by $2.02\%$ and $3.68\%$.

##  V Related Work

In the literature, prior studies have been interested in employing language
models on source code to assist software development [1, 2]. However, the
identifiers with complex names in source code make these models suffer from
the out-of-vocabulary (OOV) problem. Increasing the size of a vocabulary has a
limited effect on addressing the problem and makes models harder to scale [4].
Recently, researchers have applied open vocabulary methods for code modeling.
Karampatsis et al. [6] are the first to investigate whether open-vocabulary
methods can improve the performance of this code completion tool. They trained
a GRU-based language model with the BPE algorithm and showed that the model
performance increases over close-vocabulary models across three datasets.

At the same time, numerous works about identifier splitting have been proposed
before, and several studies have empirically compared these different
techniques [15, 16]. Some previous works also try to segment identifiers by
naming conventions [17, 18] in vocabulary construction, but no previous work
utilizes advanced identifier splitting techniques and combines them with
existing open-vocabulary algorithms.

##  VI Conclusion and Future Work

In this paper, we investigate the benefit of identifier splitting techniques
on code modeling. We propose two strategies to combine identifier splitting
with Byte-Pair Encoding (BPE) algorithm. We train open-vocabulary models with
different strategies and compare the performance over the C language dataset.
The evaluation results show that splitting identifiers improves the
performance of open-vocabulary models under a hybrid strategy, which can
improve LMs by $6.23\%$ in terms of MRR on predicting identifiers. The entropy
of LMs and R@10 on predicting identifiers are also improved by up to $2.02\%$
and $3.68\%$. Our study confirms that the potential benefits of identifier
splitting methods on open-vocabulary language models for C language.

In the future, we plan to validate our findings on more programming languages
beyond C, e.g., Java and python. Also, we are interested in considering more
models with different architectures, e.g., Transformer-based models, which
have recently drawn researchers’ attention. Besides, we plan to investigate
whether the advanced LMs can boost more code modeling-based tasks, such as
code clone detection [19], bug localization [20, 21], and code search [22].

## Acknowledgment

This research / project is supported by the National Research Foundation,
Singapore, under its Industry Alignment Fund – Pre-positioning (IAF-PP)
Funding Initiative. Any opinions, findings and conclusions or recommendations
expressed in this material are those of the authors and do not reflect the
views of National Research Foundation, Singapore.

## References

  * [1] V. Raychev, M. Vechev, and E. Yahav, “Code completion with statistical language models,” _ACM SIGPLAN Notices_ , vol. 49, 06 2014.
  * [2] Z. Chen, S. Kommrusch, M. Tufano, and et al., “Sequencer: Sequence-to-sequence learning for end-to-end program repair,” _IEEE Transactions on Software Engineering_ , vol. 47, no. 9, pp. 1943–1959, 2021.
  * [3] M. Allamanis, E. T. Barr, P. Devanbu, and C. Sutton, “A survey of machine learning for big code and naturalness,” _ACM Computing Surveys_ , vol. 51, no. 4, Jul. 2018.
  * [4] J. Xu, H. Zhou, C. Gan, Z. Zheng, and L. Li, “Vocabulary learning via optimal transport for neural machine translation,” in _2021 59th Annual Meeting of the Association for Computational Linguistics (ACL)_ , 2021, pp. 7361–7373.
  * [5] R. Sennrich, B. Haddow, and A. Birch, “Neural machine translation of rare words with subword units,” in _2016 54th Annual Meeting of the Association for Computational Linguistics (ACL)_ , 2016, pp. 1715–1725.
  * [6] R.-M. Karampatsis, H. Babii, R. Robbes, C. Sutton, and A. Janes, “Big code!= big vocabulary: Open-vocabulary models for source code,” in _2020 IEEE/ACM 42nd International Conference on Software Engineering (ICSE)_ , 2020, pp. 1073–1085.
  * [7] D. Binkley, M. Davis, D. Lawrie, and C. Morrell, “To camelcase or under_score,” in _2009 IEEE 17th International Conference on Program Comprehension (ICPC)_ , 2009, pp. 158–167.
  * [8] B. Dit, L. Guerrouj, D. Poshyvanyk, and G. Antoniol, “Can better identifier splitting techniques help feature location?” in _2011 IEEE 19th International Conference on Program Comprehension (ICPC)_ , 2011.
  * [9] E. Hill, L. Pollock, and K. Vijay-Shanker, “Automatically capturing source code context of nl-queries for software maintenance and reuse,” in _2009 IEEE 31st International Conference on Software Engineering (ICSE)_ , 2009, pp. 232–242.
  * [10] B. Vasilescu, C. Casalnuovo, and P. Devanbu, “Recovering clear, natural identifiers from obfuscated js names,” in _2017 11th Joint Meeting on Foundations of Software Engineering (ESEC/FSE)_ , 2017, p. 683–693.
  * [11] M. Hucka, “Spiral: splitters for identifiers in source code files,” _Journal of Open Source Software_ , vol. 3, p. 653, 04 2018.
  * [12] Z. Tu, Z. Su, and P. Devanbu, “On the localness of software,” in _2014 22nd ACM SIGSOFT International Symposium on Foundations of Software Engineering (FSE)_ , 2014, p. 269–280.
  * [13] K. Cho, B. van Merriënboer, C. Gulcehre, and et al., “Learning phrase representations using RNN encoder–decoder for statistical machine translation,” in _2014 19th Conference on Empirical Methods in Natural Language Processing (EMNLP)_ , 2014, pp. 1724–1734.
  * [14] V. J. Hellendoorn, S. Proksch, H. C. Gall, and A. Bacchelli, “When code completion fails: A case study on real-world completions,” in _2019 41st International Conference on Software Engineering (ICSE)_ , 2019, p. 960–970.
  * [15] E. Enslen, E. Hill, L. Pollock, and K. Vijay-Shanker, “Mining source code to automatically split identifiers for software analysis,” in _2009 IEEE 6th International Working Conference on Mining Software Repositories (MSR)_ , 2009, pp. 71–80.
  * [16] E. Hill, D. Binkley, D. Lawrie, L. Pollock, and K. Vijay-Shanker, “An empirical study of identifier splitting techniques,” _Empirical Software Engineering_ , vol. 19, 2014.
  * [17] M. Allamanis, E. T. Barr, C. Bird, and C. Sutton, “Suggesting accurate method and class names,” in _2015 10th Joint Meeting on Foundations of Software Engineering (ESEC/FSE)_ , 2015, p. 38–49.
  * [18] U. Alon, O. Levy, and E. Yahav, “code2seq: Generating sequences from structured representations of code,” in _2019 7th International Conference on Learning Representations (ICLR)_ , 2019.
  * [19] N. Mehrotra, N. Agarwal, P. Gupta, S. Anand, D. Lo, and R. Purandare, “Modeling functional similarity in source code with graph-based siamese networks,” _IEEE Transactions on Software Engineering_ , pp. 1–1, 2021.
  * [20] S. Wang and D. Lo, “Version history, similar report, and structure: Putting them together for improved bug localization,” in _2014 22nd International Conference on Program Comprehension (ICPC)_ , 2014, p. 53–63.
  * [21] Z. Yang, J. Shi, S. Wang, and D. Lo, “Incbl: Incremental bug localization,” in _2021 36th IEEE/ACM International Conference on Automated Software Engineering (ASE)_ , 2021.
  * [22] M. Lu, X. Sun, S. Wang, D. Lo, and Y. Duan, “Query expansion via wordnet for effective code search,” in _2015 IEEE 22nd International Conference on Software Analysis, Evolution, and Reengineering (SANER)_ , 2015, pp. 545–549.