# A statistical analysis of the first stages of freezing and melting of
Lennard-Jones particles: Number and size distributions of transient nuclei

Patrice Porion<EMAIL_ADDRESS>ICMN, CNRS, Université
d’Orléans, 1b rue de la Férollerie, CS 40059, 45071 Orléans cedex 02, France
Joël Puibasset<EMAIL_ADDRESS>ICMN, CNRS, Université
d’Orléans, 1b rue de la Férollerie, CS 40059, 45071 Orléans cedex 02, France

###### Abstract

The freezing/melting transition is at the heart of many natural and industrial
processes. In the classical picture, the transition proceeds via the
nucleation of the new phase, which has to overcome a barrier associated to the
free energy cost of the growing nucleus. The total nucleation rate is also
influenced by a kinetic factor which somehow depends on the number of attempts
to create a nucleus, that translates into a significant density of proto-
nuclei in the system. These transient tiny nuclei are not accessible to
experiments, but they can be observed in molecular simulations, and their
number and size distributions can be acquired and analysed. The number
distributions are carefully characterized as a function of the system size,
showing limited finite size effects. It is also shown that the proto-nuclei do
exist even in the stable phase, in agreement with the fact that the volume
contribution to their free energy is negligible in the first stages of
nucleation. Moreover, the number and size distributions vary continuously
between the stable and the metastable phases, in particular when crossing the
coexistence temperature. The size distributions associated to _any_ nucleus
and to the _largest_ one have also been calculated, and their relation
recently established for bubbles in a liquid [J. Puibasset, J. Chem. Phys.
157, 191102 (2022)] has been shown to apply here. This is an important
relation for free energy barrier calculations with biased molecular
simulations.

††preprint: AIP/123-QED

##  I Introduction

The freezing/melting transition is very common, but still not fully
understood. The main difference with the liquid/vapor transition is the fact
that it is qualitative (order/disorder transition) instead of quantitative
(high vs low fluid density). This transition is of primary importance to
understand different phenomena in nature and industry (e.g. cloud formation in
climate modeling, lubrication, tribology, cryopreservation or
reactivity)Baskaran  _et al._ (2021); Cziczo  _et al._ (2013); Lohmann and
Feichter (2005) but also regarding fundamental concepts on order/disorder
transitions (e.g. dynamics of the transition, effect of confinement or
premelting of surfaces) Slater and Michaelides (2019); Köster, Mausbach, and
Vrabec (2017); Alba-Simionesco  _et al._ (2006).

Molecular simulation is a powerful tool to study the freezing/melting
transition, in particular the molecular scale mechanisms underlying the
transition Baidakov, Bobrov, and Teterin (2011); Peng, Morris, and Aga (2010);
Trudu, Donadio, and Parrinello (2006); Gispen and Dijkstra (2023, 2024). While
in experiments melting always occurs at equilibrium by propagation from the
outer surface where premelting occurs, in simulations, the introduction of
periodic boundary conditions prevents premelting. Therefore, it is possible to
numerically monitor the experimentally inaccessible metastable (superheated)
crystal, enabling to study the underlying mechanisms of bulk crystal
destabilization during melting.Gómez  _et al._ (2003); Köster, Mausbach, and
Vrabec (2017).

Molecular simulation is also a powerful tool to study the first stages of
nucleation mechanism, in particular the activated nucleation of the new phase,
inaccessible to experiments. Measuring the equilibrium probability $p_{a}(s)$
that _any_ nucleus reaches the size $s$ gives access to the free energy
profile $W(s)$ of the embryos through Reiss, Frisch, and Lebowitz (1959);
Reiss and Bowles (1999) :

$W(s)=-kT\ln p_{a}(s)$ (1)

where $k$ is Boltzmann’s constant and $T$ is the temperature. This profile
also gives access to the free energy barrier $W(s_{c})-W(0)$ entering the
nucleation rate as given by the classical nucleation theory (CNT), where
$s_{c}$ is the size of the critical nucleus Fisher (1948); Blander and Katz
(1975); Debenedetti (1996).

Although $p_{a}(s)$ is in principle measurable in equilibrium molecular
simulations, in practice, the simulations have to be biased to sample the
regions where the free energy profile exceeds several $kT$. The bias consists
in monitoring the nucleus size from zero up to the critical size $s_{c}$
corresponding to the maximum of $W(s)$. This requires to define an order
parameter, combined with the umbrella sampling scheme of Torrie and Valleau
Torrie and Valleau (1974). In a seminal work ten Wolde, Ruiz-Montero, and
Frenkel (1996); ten Wolde and Frenkel (1997), it was proposed to chose the
size of the _largest_ nucleus as the order parameter. The corresponding
probability distribution will be denoted $p_{l}(s)$. For large nuclei, the
obtained distribution should correspond to the expected $p_{a}(s)$ since the
probability to have more than one large nucleus in the simulation box is
negligible ten Wolde, Ruiz-Montero, and Frenkel (1996); ten Wolde and Frenkel
(1997). However, for small sizes, one expects large differences between
$p_{l}(s)$ and $p_{a}(s)$. This can be seen by considering the $s\rightarrow
0$ limit: $p_{a}(s\rightarrow 0)$ is expected to become large (since the free
energy cost to produce small nuclei is negligible), while $p_{l}(s\rightarrow
0)$ is expected to become negligible (since $p_{l}(0)$ also corresponds to the
probability to have no nucleus in the metastable phase, which is highly
improbable). Furthermore, for an homogeneous fluid, $p_{a}$ is independent of
the system size, which is consistent with Eq 1, while $p_{l}$, being a maximum
of a sample, strongly depends on the total number of nuclei in the system, and
therefore depends on the system size.Casella and Berger (2002) As recognized
later, using $p_{l}$ instead of $p_{a}$ in Eq 1 may lead to
misinterpretations. Saika-Voivod, Poole, and Bowles (2006); Maibaum (2008);
Chakrabarty, Santra, and Bagchi (2008); Lundrigan and Saika-Voivod (2009);
Wedekind  _et al._ (2009); González  _et al._ (2014); Menzl  _et al._ (2016)

To circumvent the problem, it was proposed to measure directly the exact
$p_{a}(s)$ for the small values of $s$ where the disagreement is strongest,
and to complete the distribution for large $s$ values by using $p_{l}(s)$
corrected by a multiplicative constant to ensure continuity. This approach has
however some limitations, and introduces an approximation that cannot be
quantified. Goswami  _et al._ (2021) Very recently, in the context of bubble
nucleation in a metastable liquid, a general relation was proposed that
relates $p_{a}(s)$ and $p_{l}(s)$ Puibasset (2022):

$p_{a}(s)=\frac{p_{l}(s)}{\lambda_{0}\Pi_{l}(s)}$ (2)

where $\Pi_{l}(s)$ is the cumulative distribution function of $p_{l}$ and
$\lambda_{0}$ is a constant. This approach presents the advantage that it
requires only one biased simulation series to measure accurately $p_{l}(s)$
from zero up to the critical size, which is then converted straightforwardly
into $p_{a}(s)$.

The objective of this paper is to explore the robustness of this relation in
the context of the solid-liquid transition. In particular, one has to check
the hypothesis of a Poisson law for the distribution $\phi(n)$ of the number
$n$ of spontaneous nuclei in the stable and metastable phases. This hypothesis
then implies Eq. 2, where $\lambda_{0}$ is the associated Poisson parameter. A
direct calculation of $p_{a}$ and $p_{l}$ is then performed is order to
validate the relation. Furthermore, since the solid state can frequently
present different stable or metastable structures (FCC, HCP…), it is possible
to explore different transitions starting from all possible solid structures.

The paper is divided as follows: the numerical model is first presented,
including the calculation of important parameters characterizing the structure
and thermodynamics of the fluid. The distributions of the number of nuclei
$\phi(n)$ are then established for different phases, and the Fourier formalism
is used to characterize the distributions independently of the system size.
Their evolution with temperature is also established. The paper then focuses
on the nuclei size distributions $p_{a}$ and $p_{l}$, with emphasis on their
differences, and the verification of Eq. 2 that relates $p_{l}$ to $p_{a}$.
Finally, the free energy profiles are calculated to discuss the very first
stages of nucleation.

##  II Methods

_Model._ It is well-known that the interactions in solid phases are in most
cases not well approximated by pair-wise additive potentials. For example,
this is the case for metals, where many-body potentials are better
approximations. It is however assumed that the generic properties we are
interested in here are mainly controlled by thermal activation, with a free
energy profile comprising essentially volume and surface contributions. We
therefore consider a simple atomic model based on the truncated and linearly
shifted (12-6) Lennard-Jones potential (continuous forces). The energy
parameter and the kinetic diameter are denoted $\epsilon$ and $\sigma$, and
the cutoff distance is taken equal to 3$\sigma$. From now on, the quantities
are expressed in reduced units: $\sigma$ for distances, $\epsilon$ for
energies, $\epsilon/k$ for temperatures, $\epsilon/\sigma^{3}$ for pressures,
etc.

The molecular dynamics (MD) simulations are done in the isothermal-isobaric
NPT ensemble using the Large-scale Atomic/Molecular Massively Parallel
Simulator (LAMMPS)Thompson  _et al._ (2022) with Nosé/Hoover thermostat and
barostat Frenkel and Smit (2002); Hoover (1985); Nosé (1984); Andersen (1980).
The time increment is taken equal to $10^{-3}\sqrt{m\sigma^{2}/\epsilon}$ and
configurations are sampled every 1000 steps for statistical analysis.
Different system sizes are considered, ranging from $N=$ 216 to 16384 atoms.
For pure phases simulations, the initial boxes are perfectly cubic; for two-
phases simulations, the box is twice as long in one direction (juxtaposition
of two cubic boxes). During the NPT runs, the three box lengths may either be
left free to fluctuate independently or can be coupled: the results are the
same.

_Freezing, melting and coexistence temperatures._ The phase diagram for the
Lennard-Jones fluid shows a triple point slightly below $T=0.7$, with a stable
FCC solid phase.Kataoka and Yamada (2014); Travesset (2014); Schultz and Kofke
(2018); Stephan, Staubach, and Hasse (2020); Stephan  _et al._ (2019) The HCP
structure is (slightly) metastable at that temperature, and becomes stable
below $T=0.4$. On the other hand, the BCC structure is always unstable. The
exact position of coexistence lines depends on the truncation used to
calculate the interactions and the possible long range corrections, but the
phase diagram looks essentially the same.

The present work focuses on the solid-liquid transition at low pressure, i.e.
in the vicinity of the triple point (lower bound of the solid-liquid
coexistence line). These dense phases being cohesive, they remain metastable
below that point, making it possible to observe solid-liquid coexistence well
below the triple point pressure (down to possibly negative values). In our
case, the triple point pressure being quite small (around $10^{-3}$ in reduced
units for the full Lennard-Jones), we work at zero pressure, without
introducing measurable differences compared to the triple point. In
particular, the density and the structure (FCC, HCP or the radial distribution
function for the liquid) are essentially unaffected by the small change in
pressure. For instance, the relative variation of the density is only $4\times
10^{-5}$ between the triple point and zero pressure; for comparison, the
relative density variation associated with the melting of one single atom at
$T=0.62$ in a FCC crystal comprising 500 atoms is $\sim 3\times 10^{-4}$, i.e.
seven times larger (see Fig. 1).

Figure 1:  Evolution of the density of a Lennard-Jones system (500 atoms) at
zero pressure during an ascending temperature ramp from $T=0.3$ to 0.9,
followed by a descending ramp down to 0.3 (thin arrows). The large vertical
arrows correspond to the irreversible transitions between (i) the initial
perfect FCC crystal (upper blue branch) and the liquid (lowest orange branch),
and, (ii) between the liquid and the imperfectly crystallized solid
(intermediate green branch). These transitions define the melting ($T_{m}$)
and freezing ($T_{f}$) temperatures respectively. The solid-liquid equilibrium
occurs in between at $T_{sl}$ (see text).

The only noticeable difference is the fact that, strictly speaking, the liquid
and solid phases at zero pressure are both metastable with respect to the
infinitely diluted phase (i.e. the gas phase at zero pressure). However, to
simplify the discussion, the metastability with respect to the gas phase will
be ignored, and the liquid or solid phase with the lowest free energy will be
qualified as stable.

As previously mentioned, the truncation of the potential slightly influences
the phase diagram, in particular the solid-liquid coexistence line. The
coexistence temperature at zero pressure thus needs to be evaluated. We focus
on the transition between the liquid and the FCC crystal which is the most
stable structure. Starting from a FCC crystal comprising 500 atoms with
periodic boundary conditions, MD simulations with increasing temperature shows
a melting transition around $T_{m}=0.69$ (Fig. 1). Cooling back the liquid,
one observes freezing around $T_{f}=0.45$, revealing a large freezing/melting
hysteresis. Furthermore, one observes imperfect crystallization upon cooling,
with a significant number of defects in the FCC structure resulting in a lower
density and higher energy (see Fig. 1). It is emphasized that $T_{m}$ and
$T_{f}$ are only approximately positioned since they depend on the initial
conditions and the length of the MD run as well as the system size (stochastic
transitions).

This freezing/melting hysteresis defines the largest possible temperature
range where it is possible to observe both solid and liquid phases. In the
following, the calculations will be done in a slightly smaller interval,
between $T=0.46$ and 0.68, to avoid spurious effects close to the transitions.

The temperature hysteresis is split into two sub-intervals by the solid-liquid
equilibrium temperature $T_{sl}$: below $T_{sl}$, the solid is stable while
the liquid is metastable; above $T_{sl}$ it is the opposite (see Fig. 1). The
solid-liquid equilibrium has been determined by means of the two-phase
method.Kataoka and Yamada (2014) A cubic piece of FCC crystal comprising $N/2$
atoms is put in contact with a liquid box with the same number of atoms,
producing a two-phase system comprising $N$ atoms. The liquid/solid interface
is first rapidly equilibrated with the FCC atoms kept frozen. MD runs are then
performed for the $N$ atoms in the isothermal-isobaric ensemble at $P$ = 0 at
different temperatures. In each case, the evolution of the total number of FCC
atoms is monitored. At low $T$, the liquid region invariably fully
crystallizes, while at high $T$ the crystal always melts. For intermediate
temperatures, the system may either crystallize or melt depending on the
initial conditions. One can define a probability towards crystallization by
monitoring different initial conditions and counting the number of
configurations that lead to a crystal. Figure 2 shows the evolution of this
probability with the temperature $T$ for three system sizes.

Figure 2:  Symbols: crystallization probability for a system initially
comprising an FCC crystal in contact with an equal number of liquid atoms, as
a function of the temperature, at zero pressure. Three system sizes are
considered (total number of crystal+liquid atoms $N$= 512, 1000 and 4096). The
crystallization probability is calculated based on a large number of initial
conditions: 2000 runs for $N$ = 512 and 1000; 500 runs for $N$ = 4096. Lines:
sigmoid fits. Inset: lower and upper bound for the coexistence temperature as
given by a crystallization probability between 5 and 95%.

As can be seen, the crystallization probability decreases from one, at low
temperature, down to zero, at high $T$, following a sigmoid curve. As
expected, the larger the system, the sharper the transition. The solid-liquid
equilibrium temperature $T_{sl}$ can be defined as the temperature where the
crystallization probability equals 0.5 (crystallization and melting are
equiprobable), and an interval can be defined by imposing a crystallization
probability between 5 and 95% (see inset). As can be seen, $T_{sl}$ slightly
varies with the system size, with values around 0.62, in agreement with
previous determinations based on other techniques.Baidakov and Tipeev (2012);
Malek, Morrow, and Saika-Voivod (2015); Morris and Song (2003)

_Order parameter._ The local lattice structure around an atom is determined
using the Polyhedral Template Matching (PTM) method. Larsen, Schmidt, and
Schiøtz (2016) For each atom of the sample, its local structure is compared to
that of a list of templates (FCC, HCP, BCC, etc.), where “local” means that
one includes the neighbors within the first minimum of the radial distribution
function at temperature $T$. For each atom and template, the matching level is
quantified by a scalar (denoted RMSD for root-mean-square deviation). Each
atom is then assigned the lowest RMSD value among all templates, and
temporarily tagged with the template that matches the best with its local
structure, i.e. the one giving the lowest RMSD. By construction, the RMSD is a
scalar ranging from zero to infinity. It takes well defined values for atoms
in perfect crystals (FCC, HCP…), but values belonging to enlarged
distributions when thermal agitation comes into play.

Figure 3 gives the RMSD distributions calculated for the three phases of
interest for this work (an FCC crystal with 500 atoms, an HCP crystal with 512
atoms, and a liquid sample with 500 atoms) at $P$ = 0 and $T$ = 0.62,
corresponding to the FCC-liquid coexistence temperature previously determined.

Figure 3:  RMSD distributions calculated for the three samples corresponding
to the phases of interest for this work: a FCC crystal (500 atoms), an HCP
crystal (512 atoms) and a liquid sample (500 atoms), at $P$ = 0 and $T$ =
0.62. The optimal value of the RMSD cutoff, at the intersection between the
liquid and the two solid curves, is given by the vertical dashed line at RMSD0
= 0.15.

The RMSD values for the two crystalline phases (FCC and HCP) are identically
distributed around the nominal value 0.09 with similar widths, each phase
matching best with its corresponding template (FCC or HCP), as expected. On
the other hand, the RMSD density for the liquid phase is distributed around a
larger value (0.24) with a larger spreading, and the corresponding RMSD values
correspond to FCC, HCP or BCC templates, in proportion corresponding
approximately to 10, 70 and 20% respectively.

As can be seen, the RMSD distributions for the solid phases (FCC or HCP) and
for the liquid phase are well separated. The intersection, denoted RMSD0, is
used as a threshold to identify the local structure of the atoms in any given
molecular configuration, as follows. Each atom is assigned the lowest RMSD
value among all templates. If its RMSD falls below RMSD0, the atom is tagged
as FCC or HCP, according to the corresponding RMSD optimal matching. On the
other hand, if the RMSD is larger than RMSD0, the local structure is
identified as disordered and the atom is tagged as liquid. In agreement with
other authors, we found RMSD0 = 0.15 (see Fig. 3). Chung  _et al._ (2022);
Çamkıran, Parsch, and Hibbard (2022) The PTM method is implemented in the
Ovito software used to analyse the configurations.Stukowski (2010)

_Nuclei: identification and statistics._ Let us first consider nucleation in
the liquid phase. Visual inspection shows that solid-like atoms of different
structures (FCC, HCP or BCC) appear spontaneously in the system [see Fig.
4(a)].

Figure 4:  Snapshots (generated with OvitoStukowski (2010)) illustrating the
procedure to detect the solid nuclei in the liquid (a), the defective (non-
FCC) nuclei in the FCC crystal (b), and the pure liquid nuclei in the FCC
crystal (c). Left column: all atoms are shown; right column: only the atoms
corresponding to the selected structures are shown, the nuclei being defined
as clusters of these atoms. The color code corresponds to the local structure
as given by the PTM method.

Quite remarkably, the relative abundance is in favor of HCP, at odds with the
fact that FCC is the most stable structure for the bulk solid (at the
coexistence temperature $T$ = 0.62, the HCP-to-FCC ratio equals 3.7). This
confirms that the relative stability of tiny crystallites in liquid does not
necessarily follow the bulk behavior Bulutoglu  _et al._ (2023); Desgranges
and Delhommelle (2007); Leyssale, Delhommelle, and Millot (2005). Furthermore,
these solid-like atoms, which are mostly isolated, sometimes group into
clusters that are not necessarily pure FCC or HCP, as can be seen on the
example shown in Fig. 4(a). We are therefore led to define the clusters
without specifying the particular solid structure. The nuclei are thus defined
as clusters of any solid atoms (FCC, HCP and BCC), the cluster grouping
criterion being chosen to be the largest value, among the different crystal
structures, of the first minimum of the radial distribution function at each
temperature $T$.

Let us now consider the tiny nuclei appearing in the FCC crystal. Visual
inspection shows that “liquid atoms” is not the unique kind of defects that
appear: HCP or BCC atoms are also observed [see Fig. 4(b)]. As shown,
different defective (non-FCC) atoms can appear in the same cluster: we are
therefore led to define the nuclei in the FCC crystal as clusters of non-FCC
atoms (including liquid, HCP and BCC), by analogy with what we did for the
solid nuclei in the liquid phase. Note that in this case, the distance
criterion for clustering corresponds to the first minimum of the radial
distribution function of the liquid (larger than that of the solid).

However, it is also possible to restrict the cluster analysis to liquid-like
atoms: this is what is depicted in Fig. 4(c). In that case, only liquid-like
atoms are taken into account to define the nuclei, with the corresponding
cluster criterion. These clusters will be denoted “liquid nuclei”. It is not
clear which definition of nuclei (liquid or non-FCC) is the best to describe
nucleation in solids. Obviously, one expects measurable differences for small
clusters comprising few atoms. For instance, Fig. 4(b) and 4(c) illustrate a
case where the number and the size of nuclei differ depending on the
criterion. However, for large nuclei relevant to calculate nucleation
barriers, their composition is expected to be dominated by the most stable
phase, and therefore, both criteria should, at the end, give similar results.

For the HCP crystal, the situation is very similar to the FCC case: the nuclei
may either be defined as clusters of non-HCP atoms (liquid, FCC and BCC), or
clusters of liquid atoms solely.

This procedure allows to define the nuclei appearing in the system for each
atomic configuration, and determine their number $n$ and their sizes $s$
defined as the total number of atoms in each cluster. The distribution of the
number of nuclei per configuration is denoted $\phi(n)$; it is normalized to
one. The densities $p_{a}(s)$ and $p_{l}(s)$ resume to discrete histograms
since the nucleus size $s$ is defined as its number of atoms. Note that each
configuration contributes once for $p_{l}$, leading to a histogram normalized
to one. On the other hand, each configuration contributes as many times to
$p_{a}$ as the number of nuclei in the configuration. Therefore, after
building the histogram over all configurations, it has to be divided by the
number of configurations times the average number of nuclei in the system, so
that it is normalized to one. The number of configurations used for the
averages ranges from 20000 for the largest systems to $10^{6}$ for the
smallest ones. Note also that $p_{a}(s)$ is defined only for non-zero sizes,
while $p_{l}(s)$ is defined for any size, including zero; in particular,
$p_{l}(0)$ is the probability that the molecular configuration contains no
nuclei.

##  III Results

###  III.1 Distribution of the number of nuclei $\phi(n)$

_Influence of the system size._ Typical $\phi_{N}(n)$ (in this section and the
next, the $N$-dependence is shown as a subscript) are given in Fig. 5 in the
case of solid nucleation in the liquid, at the solid-liquid coexistence
temperature $T=0.62$ and for different system sizes from $N$ = 250 to 16000
atoms.

Figure 5:  Symbols: distributions $\phi_{N}(n)$ of the number $n$ of solid
nuclei (FCC+HCP+BCC) in the liquid at $P=0$ and $T=0.62$, shown in log-linear
scale and for different system sizes $N$ given in the figure. Lines are fits
with Poisson distributions.

The distributions exhibit a maximum, corresponding to a typical number of
nuclei, that increases linearly with the system size. The width of the
distributions also increases with the system size. One expects that nucleation
should follow a Poisson distribution, at least in the low rate regime where
the nucleation events are spatially independent.Maibaum (2008); Blow  _et al._
(2023) This is confirmed by the good quality of the Poisson fits shown in Fig.
5.

The symmetric situation, corresponding to the nucleation in the FCC crystal at
the same temperature and pressure, gives similar results shown in Fig. 6.

Figure 6:  Symbols: distributions $\phi_{N}(n)$ of the number $n$ of nuclei in
the FCC crystal at $P=0$ and $T=0.62$, shown in log-linear scale and for
different system sizes $N$ given in the figure. The empty symbols refer to the
non-FCC nuclei, while the stars refer to the pure liquid nuclei (shown only
for the largest system). Lines are fits of the non-FCC nuclei distributions
with the Poisson law.

For comparison, the distributions for the largest system are calculated either
with all non-FCC (liquid+HCP+BCC) atoms to define the clusters, or with only
liquid atoms. The main difference appears as a shift of the distribution to
slightly lower numbers of nuclei for pure liquid clusters (see Fig. 6).
However, the relative difference is less than 3%. Furthermore, this difference
is expected to decrease for larger nuclei (closer to the nucleation barrier).
Focusing on non-FCC nuclei, the behavior with the system size is essentially
the same as for solid nuclei in liquid, the main difference being the larger
values for the average number of nuclei, showing that, at the same temperature
and pressure, more nuclei are detected in the solid compared to the liquid.
The second important difference is that the Poisson fits exhibit small
disagreements, mainly for the lowest values of the distributions (see Fig. 6).

The same system size analysis has been done for the HCP crystal, showing
distributions almost indistinguishable from those obtained in the FCC crystal
for similar system sizes (for instance 500 and 512 or 4000 and 4096 in the FCC
and HCP crystals respectively). We observe in particular the same departure
from the Poisson law. In order to elucidate whether this non-Poisson behavior
is due to finite size effects, we need to compare quantitatively the
distributions associated to different system sizes, as follows.

_Finite size analysis._ To characterize the evolution of the distributions
with the system size, it is useful to notice that the distribution for the
system of size $N_{1}+N_{2}$ is the convolution of the distributions for the
sub-systems of size $N_{1}$ and $N_{2}$. For instance, the probability to
observe no nuclei in the total system is
$\phi_{N_{1}+N_{2}}(0)=\phi_{N_{1}}(0)\phi_{N_{2}}(0)$. More generally, one
has $\phi_{N_{1}+N_{2}}(n)=\sum_{i=0}^{n}\phi_{N_{1}}(i)\phi_{N_{2}}(n-i)$.

Therefore, introducing the Fourier formalism for $\phi_{N}(n)$:

$\tilde{\phi}_{N}(q)=\sum_{n=0}^{\infty}\phi_{N}(n)e^{iqn}$ (3)

and

$\phi_{N}(n)=\frac{1}{2\pi}\int_{-\pi}^{\pi}\tilde{\phi}_{N}(q)e^{-iqn}dq$ (4)

one has the remarkable property:

$\tilde{\phi}_{N_{1}+N_{2}}(q)=\tilde{\phi}_{N_{1}}(q)\tilde{\phi}_{N_{2}}(q)$
(5)

meaning that $\ln\tilde{\phi}_{N}(q)$ is extensive with the system size $N$.
We therefore introduce the normalized quantity

$\psi_{N}(q)=N^{-1}\ln\tilde{\phi}_{N}(q)$ (6)

to characterize and compare the distributions for different system sizes. This
quantity is shown in Fig. 7 in the case of nucleation in the FCC crystal at
$P=0$ and $T=0.62$.

Figure 7:  Real and imaginary parts of $\psi_{N}(q)$ (Eq 6) for the nuclei
distributions $\phi_{N}(n)$ (shown in Fig. 6) in the FCC crystal at $P=0$ and
$T=0.62$ and for different system sizes given in the figure. Upper panel: real
part versus $q$ (main) and $q^{2}$ (insert). Lower panel: imaginary part
versus $q$. The small insets are local enlargements to evidence the system
size dependence.

The results for nucleation in HCP are essentially indistinguishable from FCC,
and, for nucleation in liquid, they are qualitatively the same except that the
system size dependence is not visible. For these reasons they are not shown.
As evidenced in Fig. 7, the imaginary part is essentially linear with the
frequency $q$, while the real part is quadratic (see inset). A best fit with
the equation

$\psi_{N}(q)=ia_{N}q-0.5b_{N}q^{2}$ (7)

gives the coefficients $a_{N}$ and $b_{N}$ (corresponding to the first
cumulants) as a function of the system size $N$. As shown in Fig. 8, $a_{N}$
and $b_{N}$ vary linearly with $1/N$, as expected in standard finite size
analysis, suggesting the values $a_{\infty}=9.23\times 10^{-3}$ and
$b_{\infty}=8.02\times 10^{-3}$ for nucleation in the infinite FCC crystal
(quantitatively the same result is obtained for the HCP structure).

Figure 8:  Linear ($a_{N}$, down triangles) and quadratic ($b_{N}$, up
triangles) coefficients (Eq 7) for, respectively, the imaginary and real parts
of $\psi_{N}(q)$ (see Fig. 7), as a function of $1/N$, for nucleation in the
liquid (blue points) and in the FCC (green) and HCP (red) crystals. A best fit
is shown with one (liquid) or two parallel (FCC crystal) lines. Circles:
$\lambda_{N}=-N^{-1}\ln\phi_{N}(0)$ when defined with a good accuracy (i.e.
for the smallest systems).

On the other hand, the coefficient $a_{N}$ and $b_{N}$ associated to
nucleation in the liquid phase exhibit no measurable dependence with $N$ and
are almost identical: $a_{\infty}=3.22\times 10^{-3}$ and
$b_{\infty}=3.20\times 10^{-3}$.

For independent nucleations, one expects that $\phi_{N}$ should follow a
Poisson distribution $\phi^{P}_{N}(n)=e^{-\Lambda_{N}}\Lambda_{N}^{n}/n!$,
with a parameter $\Lambda_{N}$ (equal to the average number of nuclei in the
system) proportional to the system size: $\Lambda_{N}=N\lambda$. In this case,

$\displaystyle\psi_{N}(q)=N^{-1}\ln\tilde{\phi}^{P}_{N}(q)=\lambda\left(e^{iq}-1\right)$
$\displaystyle=i\lambda q-0.5\lambda q^{2}+O(q^{3})$ (8)

meaning that the coefficients $a_{N}=b_{N}=\lambda$. This is what we observe
for nucleation in the liquid phase whatever the size of the system (see Fig.
8). The fact that in the FCC and the HCP crystals the two coefficients $a_{N}$
and $b_{N}$ slightly differ (by a constant independent of the system size as
shown by the two parallel lines in Fig. 8) suggests that the distributions of
the number of nuclei slightly depart from the Poisson law even for the
infinite system. The disagreement is however moderate ($<10\%$ in the lowest
part of the distributions where the disagreement is strongest, see Fig. 6). It
is not yet understood why the distributions depart from the Poisson
statistics, in a way independent of the system size. It is suggested that
correlations between nuclei could explain this discrepancy.

Note that for the Poisson law $\phi^{P}_{N}$,
$\lambda=-N^{-1}\ln\phi^{P}_{N}(0)$. By analogy, for any $\phi_{N}$, we denote
$\lambda_{N}=-N^{-1}\ln\phi_{N}(0)$, which now depends on the system size $N$.
Figure 8 shows that $\lambda_{N}$ (circles) is an excellent guess for the
coefficient $a_{\infty}$ of the infinite system, at least as long as the
system is small enough so that $\phi_{N}(0)$ can be determined with a good
accuracy. It is however not yet understood why this quantity, that can be
determined only for the smallest systems, is a good guess for $a_{\infty}$.

This finite size analysis shows that $N\geq 2000$ is a good compromise between
computational cost and accuracy for all systems (less than $1\%$ discrepancy
compared to the infinite system). In the following, the system size is imposed
and explicitly given for each situation: the subscripts $N$ will be dropped.

_Influence of the temperature and the nature of the transition._ The
distributions $\phi(n)$ of the number of solid nuclei in the liquid are shown
in Fig. 9 for three temperatures (close to $T_{f}$, $T_{m}$ and in between).

Figure 9:  Symbols: distributions $\phi(n)$ of the number of solid
(FCC+HCP+BCC) nuclei in the liquid system with $N=2000$ at $P=0$ and for
different temperatures given in the figure. The lines are fits with Poisson
distributions.

As can be seen, the lower the temperature, the larger the probability to
observe a large number of nuclei and the lower the probability to observe a
small number of nuclei. It is also visible that the Poisson fits are less
accurate at low temperature, suggesting correlations between nuclei when their
density increases.

The inverse situation has been studied, corresponding to nucleation in the FCC
crystal. The distributions of the number of non-FCC nuclei are given in Fig.
10 for different temperatures.

Figure 10:  Symbols: distributions $\phi(n)$ of the number of non-FCC
(liquid+HCP+BCC) nuclei in the FCC crystal of size $N=2048$ at $P=0$ and for
different temperatures given in the figure. The lines are fits with Poisson
distributions.

The trend with temperature is now reversed, with the average number of nuclei
increasing for larger $T$. The Poisson fits also exhibit good matching at low
temperature, while the highest temperature clearly exhibits departure from the
expected Poisson distribution for independent events. As previously, this
correlates with the increase in the number of nucleation events.
Quantitatively similar results are obtained for nucleation in the HCP crystal.

As previously, these distributions have been characterized through the two
coefficients $a$ and $b$ giving the linear and quadratic dependence of
$\psi(q)$ (Eq. 6 and 7). Figure 11 gives their evolution with temperature.

Figure 11:  Evolution with $T$ of the coefficients $a$ and $b$ giving the
linear and quadratic terms of the development of $\psi(q)$ at small $q$ (Eq. 6
and 7). From top to bottom: nucleation in the liquid ($N=2000$), in the FCC
crystal ($N=2048$) and in the HCP crystal ($N=2744$) respectively. In each
case, several criteria have been considered to define the nuclei (given in the
figures). The hatched region in the liquid and FCC panels corresponds to the
temperature range where coexistence between liquid and FCC is likely to occur
for a system of size $N\simeq 2000$ (see Fig 2).

The first point to be noticed is that $a$ and $b$ decrease when temperature
increases for solid nucleation in liquid, while it it the converse for
nucleation in the two solids. It is also visible that the crystal structure
(FCC vs HCP) has a negligible impact on $a$ and $b$ (absolute differences less
than $10^{-4}$). Therefore, the crystal structure is shown to have a
negligible influence on the first stages of nucleation.

Focusing on the FCC crystal, the data show that the two definitions of the
nuclei (pure liquid vs non-FCC) give essentially the same results. One can
however notice few tiny differences, not imputable to statistical fluctuations
since very similar features are observed for the HCP crystal. For nucleation
in the liquid, in addition to the usual definition including all solids, the
distributions associated with pure FCC, HCP or BCC nuclei were also analysed.
The main point to be noticed is that the BCC nuclei are the less abundant,
while HCP contribution is the largest, at odds with the fact that the FCC
phase is the most stable.

Focusing on solid nucleation in liquid, one finds $a=b$ for $T\geq 0.58$,
while this equality fails below 0.56, the discrepancy increasing at lower
temperatures. This means that the distribution $\phi(n)$ does not follow the
Poisson law at low $T$, revealing a lack of independence between nucleation
events. Examination of nucleation in the two solids (FCC and HCP) shows a
symmetric situation where $a=b$ below $T=0.58$ and departure from the Poisson
law above $T=0.59$. Note that, in this case, the difference between $a$ and
$b$ grows faster than for solid nucleation in the liquid. As is clearly
evidenced by the point at $T=0.61$, the difference between $a$ and $b$ appears
strictly below the coexistence temperature $T_{sl}$ and is therefore not
related to it. Furthermore, one does not observe any discontinuity or change
in slope while crossing the solid-liquid coexistence temperature: the change
in stability of the phase (stable vs metastable) has no impact on the
thermodynamic properties relevant for the first stages of nucleation. As a
matter of fact, the first stages are dominated by surface contributions, not
influenced by the stability of the phase. On the other hand, for larger
nuclei, where volume contributions play a role, the relative stability between
the phases is expected to become important.

It is also interesting to examine the data close to the freezing and melting
temperatures. One expects divergence of the number of nuclei while approaching
the spinodals (where nucleation barriers are expected to reach zero). The
coefficient $a$, giving the average number of nuclei, is thus expected to
increase faster for the liquid around $T_{f}$ or for the solids around
$T_{m}$. This is qualitatively observed for solid nucleation in the liquid. On
the other hand, nucleation in the two crystals exhibit an inflexion point
below $T_{m}$, suggesting a saturation effect. Exploring higher temperatures,
however, is made more difficult due to the rapid melting of the system.

###  III.2 Distributions of nucleus size $p_{a}(s)$ and $p_{l}(s)$

_Any vs largest nuclei._ The histograms giving the size distributions of any
$p_{a}(s)$ or the largest $p_{l}(s)$ solid nuclei in the liquid at $T$ = 0.62
are given in Fig. 12 for different system sizes ranging from $N$ = 250 to
4000.

Figure 12:  Symbols: size distributions of any ($p_{a}(s)$, upper panel) and
of the largest ($p_{l}(s)$, lower panel) solid (FCC+HCP+BCC) nuclei in the
liquid at $P=0$ and $T=0.62$, for different system sizes $N$ given in the
figure. Dotted lines (lower panel) are guides to the eye. Lines (upper panel):
calculation of $p_{a}(s)$ from $p_{l}(s)$, as given by Eq 9, following the
same color code. Inset: $-\ln p_{l}(0)$ vs $N$.

As can be seen, $p_{a}(s)$ is independent of the system size, since all nuclei
contribute equally to the distribution, independently of the others. It is an
intrinsic statistical property of the fluid, that depends only on intensive
parameters like $T$ or $P$, relevant for nucleation processes, as given by Eq
1.

On the other hand, $p_{l}(s)$ is system size dependent. The larger the system,
the lower the probability to observe a largest nucleus smaller than few atoms.
The most stringent situation is for $p_{l}(0)$ which corresponds to the
probability not to observe nuclei in the system: it decreases exponentially
with the system size (see inset). Conversely, the probability to observe a
largest nucleus with a large number of atoms increases with the system size.

The results of nucleation in the FCC crystal at the same temperature $T$ =
0.62 are given in Fig. 13.

Figure 13:  Symbols: size distributions of any ($p_{a}(s)$, upper panel) and
of the largest ($p_{l}(s)$, lower panel) nuclei in the FCC and HCP crystals at
$P=0$ and $T=0.62$, for different system sizes $N$ given in the figure. Dotted
lines (lower panel) are guides to the eye. Lines (upper panel): calculation of
$p_{a}(s)$ from $p_{l}(s)$, as given by Eq 9, following the same color code.
Inset: $-\ln p_{l}(0)$ vs $N$.

As can be seen, the same qualitative behavior is observed for the FCC crystal,
i.e. a large system size dependence for $p_{l}$ while $p_{a}$ remains
unchanged. The only difference is quantitative: the probability to observe a
given nucleus size is significantly lower for solid in liquid compared to the
converse at the same coexistence temperature $T=0.62$.

The HCP crystal exhibits quantitatively similar results to the FCC case. This
is clear for $p_{a}$ from Fig. 13. For $p_{l}$, since it depends on the system
size, a direct comparison between FCC and HCP requires to choose the system
sizes closest as possible: this is done for $N_{\rm HCP}=512$ which exhibits
quantitatively the same distributions as for $N_{\rm FCC}=500$ (see Fig. 13).

_Conversion of $p_{l}$ into $p_{a}$._ As mentioned in the Introduction,
$p_{a}(s)$ is not accessible to simulations for large $s$, while $p_{l}(s)$ is
calculable thanks to biased simulations for any size $s$. In Ref. [Puibasset,
2022] (Eq. 12) it was proposed the following algorithm to convert the discrete
$p_{l}(s)$ into $p_{a}(s)$ which applies in our case since the distributions
$\phi$ are close to Poisson law:

$p_{a}(s)=\frac{1}{\lambda}\left\\{\ln\Pi_{l}(s)-\ln\Pi_{l}(s-1)\right\\}\text{
for }s\geq 1$ (9)

where $\Pi_{l}(s)=\sum_{i=0}^{s}p_{l}(i)$ is the discrete cumulative
distribution of $p_{l}$ and $\lambda=-\ln p_{l}(0)$. The results are given in
the upper panels of Figs. 12 and 13 as solid lines. As can be seen, the
results perfectly superimpose to the data $p_{a}(s)$ established by direct
simulations, which validates the transformation.

##  IV Discussion and conclusion.

The objective of this work was to validate Eq 9 in the case of the solid-
liquid transition. In order to discuss nucleation, it is now required to
calculate $p_{l}(s)$ for any nucleus size, in particular up to the top of the
hill corresponding to the nucleation barrier. This requires lengthy biased
simulations left for a further work. It is however already possible to discuss
the very first stages of nucleation. In particular, using Eq 1, the free
energy profiles $W(s)$ can be calculated (up to an irrelevant constant): they
are shown in Figs. 14 and 15 for nucleation in the liquid and FCC crystal
respectively.

Figure 14:  Main panel: Symbols: free energy profiles of solid nucleation in
the liquid as given by Eq 1 for different temperatures and for $N=2000$. Lines
are best fits with Eq 10; the evolution of the free parameters $\gamma$ and
$\delta$ with $T$ is given in the insets.  Figure 15:  Main panel: Symbols:
free energy profiles of nucleation in the FCC crystal as given by Eq 1 for
different temperatures and for $N=2048$. Lines are best fits with Eq 10; the
evolution of the free parameters $\gamma$ and $\delta$ with $T$ is given in
the insets.

As can be seen, the free energy profiles are monotonously increasing with the
nucleus size. Their evolution with temperature is such that the free energy
profile is larger in the stable phase (liquid at high $T$ or FCC crystal at
low $T$) and lower in the metastable phase (liquid at low $T$ or FCC crystal
at high $T$). In particular, for the most metastable cases, the free energy
profiles suggest the existence of a maximum that was however out of accessible
data.

In the classical nucleation theory, the free energy profile is that of a
growing nucleus with a perfectly spherical shape. The main contributions are
homogeneous to surface and volume, corresponding respectively to the free
energy of the interface and that of the difference between the old and the new
phases. We thus introduce the general equation:

$W(s)-W(0)=\gamma\hskip 1.0pts^{2/3}+\delta s\hskip 1.0pt$ (10)

where the free parameter $\gamma$ plays the role of an effective interface
free energy and $\delta$ plays the role of a free energy difference between
the phases. A best fit of the data with Eq 10 is shown in Figs. 14 and 15 as
lines, and the evolution of the free parameters with temperature is given in
the insets. The first point to be mentioned is that the parameters vary
continuously with temperature, showing no discontinuity while crossing the
coexistence temperature ($T$ = 0.62). As expected, the surface parameter is
always positive. However, it takes different values depending on whether we
are considering solid nucleation in the liquid or the converse, and it
increases with increasing temperature. For nucleation in the FCC crystal, the
$\delta$ parameter roughly decreases from a positive value below the
coexistence temperature to a negative value above that temperature: this
behavior is in agreement with what is expected from a free energy difference
between the bulk phases. However, for solid nucleation in liquid, the behavior
is at odds with expectations: the free energy difference is always negative.
It is however emphasized that we focused on the very first stages of
nucleation where the volume contribution is small compared to the surface
therm. Further calculations (at larger nucleus size) are required to discuss
with a better accuracy the evolution of the parameter $\delta$ with $T$.

As a conclusion, the careful analysis of the number distributions of nuclei
$\phi(n)$ in the stable and metastable phases exhibit no qualitative
differences. The distributions follow essentially the expected Poisson law for
independent nucleations, except in the case where the density of nuclei is so
large that correlations may appear. This is in agreement with the results of
the finite size analysis showing no size dependence at low nuclei density. In
this context, the relation (Eq 9) between the size distributions of any
($p_{a}$) or the largest ($p_{l}$) nucleus is expected to hold, and this is
what has been shown in the liquid as well as in the FCC and HCP crystal
phases. The main objective is now to use biased molecular simulations to
establish $p_{l}$ on the whole range of nucleus size in order to discus the
evolution of the nucleation profiles with the temperature.

###### Acknowledgements.

The authors are grateful to the CNRS interdisciplinary “Défi Needs” through
its “MiPor” program (Project DARIUS), and to the financial support of Agence
Nationale de la Recherche through the projects ANR-17-CE30-0002 (CavConf) and
ANR-23-CE30-0028-03 (NanoCav).

## Data Availability Statement

The data that support the findings of this study are available from the
corresponding author upon reasonable request.

## References

## References

  * Baskaran  _et al._ (2021) A. Baskaran, M. Kaari, G. Venugopal, R. Manikkam, J. Joseph, and P. V. Bhaskar, “Anti freeze proteins (Afp): Properties, sources and applications – A review,” Int. J. Biol. Macromol. 189, 292–305 (2021).
  * Cziczo  _et al._ (2013) D. J. Cziczo, K. D. Froyd, C. Hoose, E. J. Jensen, M. Diao, M. A. Zondlo, J. B. Smith, C. H. Twohy, and D. M. Murphy, “Clarifying the dominant sources and mechanisms of cirrus cloud formation,” Science 340, 1320–1324 (2013).
  * Lohmann and Feichter (2005) U. Lohmann and J. Feichter, “Global indirect aerosol effects: a review,” Atmos. Chem. Phys. 5, 715–737 (2005).
  * Slater and Michaelides (2019) B. Slater and A. Michaelides, “Surface premelting of water ice,” Nat. Rev. Chem. 3, 172–188 (2019).
  * Köster, Mausbach, and Vrabec (2017) A. Köster, P. Mausbach, and J. Vrabec, “Premelting, solid-fluid equilibria, and thermodynamic properties in the high density region based on the Lennard-Jones potential,” J. Chem. Phys. 147, 144502 (2017).
  * Alba-Simionesco  _et al._ (2006) C. Alba-Simionesco, B. Coasne, G. Dosseh, G. Dudziak, K. E. Gubbins, R. Radhakrishnan, and M. Sliwinska-Bartkowiak, “Effects of confinement on freezing and melting,” J Phys: Condens Matter. 18, R15–R68 (2006).
  * Baidakov, Bobrov, and Teterin (2011) V. G. Baidakov, K. Bobrov, and A. S. Teterin, “Cavitation and crystallization in a metastable Lennard-Jones liquid at negative pressures and low temperatures,” J. Chem. Phys. 135, 054512 (2011).
  * Peng, Morris, and Aga (2010) L. J. Peng, J. R. Morris, and R. S. Aga, “A parameter-free prediction of simulated crystal nucleation times in the Lennard-Jones system: From the steady-state nucleation to the transient time regime,” J. Chem. Phys. 133, 084505 (2010).
  * Trudu, Donadio, and Parrinello (2006) F. Trudu, D. Donadio, and M. Parrinello, “Freezing of a Lennard-Jones fluid: From nucleation to spinodal regime,” Phys. Rev. Lett. 97, 105701 (2006).
  * Gispen and Dijkstra (2023) W. Gispen and M. Dijkstra, “Brute-force nucleation rates of hard spheres compared with rare-event methods and classical nucleation theory,” J. Chem. Phys. 159, 086101 (2023).
  * Gispen and Dijkstra (2024) W. Gispen and M. Dijkstra, “Finding the differences: Classical nucleation perspective on homogeneous melting and freezing of hard spheres,” J. Chem. Phys. 160, 141102 (2024).
  * Gómez  _et al._ (2003) L. Gómez, A. Dobry, C. Geuting, H. T. Diep, and L. Burakovsky, “Dislocation lines as the precursor of the melting of crystalline solids observed in monte carlo simulations,” Phys. Rev. Lett. 90, 095701 (2003).
  * Reiss, Frisch, and Lebowitz (1959) H. Reiss, H. L. Frisch, and J. L. Lebowitz, “Statistical mechanics of rigid spheres,” J. Chem. Phys. 31, 369–380 (1959).
  * Reiss and Bowles (1999) H. Reiss and R. K. Bowles, “Some fundamental statistical mechanical relations concerning physical clusters of interest to nucleation theory,” J. Chem. Phys. 111, 7501–7504 (1999).
  * Fisher (1948) J. C. Fisher, “The fracture of liquids,” J. Appl. Phys. 19, 1062–1067 (1948).
  * Blander and Katz (1975) M. Blander and J. L. Katz, “Bubble nucleation in liquids,” AIChE Journal 21, 833–848 (1975).
  * Debenedetti (1996) P. G. Debenedetti,  _Metastable Liquids: Concepts and Principles_ (Princeton University Press, Princeton, NJ, 1996).
  * Torrie and Valleau (1974) G. M. Torrie and J. P. Valleau, “Monte carlo free energy estimates using non-Boltzmann sampling: Application to the sub-critical Lennard-Jones fluid,” Chem. Phys. Lett. 28, 578–581 (1974).
  * ten Wolde, Ruiz-Montero, and Frenkel (1996) P.-R. ten Wolde, M. J. Ruiz-Montero, and D. Frenkel, “Simulation of homogeneous crystal nucleation close to coexistence,” Faraday Discussions 104, 93–110 (1996).
  * ten Wolde and Frenkel (1997) P. R. ten Wolde and D. Frenkel, “Enhancement of protein crystal nucleation by critical density fluctuations,” Science 277, 1975–1978 (1997).
  * Casella and Berger (2002) G. Casella and R. L. Berger,  _Statistical Inference_ (Thomson Learning, New York, 2002).
  * Saika-Voivod, Poole, and Bowles (2006) I. Saika-Voivod, P. H. Poole, and R. K. Bowles, “Test of classical nucleation theory on deeply supercooled high-pressure simulated silica,” J. Chem. Phys. 124, 224709 (2006).
  * Maibaum (2008) L. Maibaum, “Comment on “Elucidating the mechanism of nucleation near the gas-liquid spinodal”,” Phys. Rev. Lett. 101, 019601 (2008).
  * Chakrabarty, Santra, and Bagchi (2008) S. Chakrabarty, M. Santra, and B. Bagchi, “Chakrabarty, Santra, and Bagchi reply,” Phys. Rev. Lett. 101, 019602 (2008).
  * Lundrigan and Saika-Voivod (2009) S. E. M. Lundrigan and I. Saika-Voivod, “Test of classical nucleation theory and mean first-passage time formalism on crystallization in the Lennard-Jones liquid,” J. Chem. Phys. 131, 104503 (2009).
  * Wedekind  _et al._ (2009) J. Wedekind, G. Chkonia, J. Wölk, R. Strey, and D. Reguera, “Crossover from nucleation to spinodal decomposition in a condensing vapor,” J. Chem. Phys. 131, 114506 (2009).
  * González  _et al._ (2014) M. A. González, G. Menzl, J. L. Aragones, P. Geiger, F. Caupin, J. L. F. Abascal, C. Dellago, and C. Valeriani, “Detecting vapour bubbles in simulations of metastable water,” J. Chem. Phys. 141, 18C511 (2014).
  * Menzl  _et al._ (2016) G. Menzl, M. A. Gonzalez, P. Geiger, F. Caupin, J. L. F. Abascal, C. Valeriani, and C. Dellago, “Molecular mechanism for cavitation in water under tension,” PNAS 113, 13582 (2016).
  * Goswami  _et al._ (2021) Y. Goswami, V. V. Vasisht, D. Frenkel, P. G. Debenedetti, and S. Sastry, “Thermodynamics and kinetics of crystallization in deeply supercooled Stillinger–Weber silicon,” J. Chem. Phys. 155, 194502 (2021).
  * Puibasset (2022) J. Puibasset, “A general relation between the largest nucleus and all nuclei distributions for free energy calculations,” J. Chem. Phys. 157, 191102 (2022).
  * Thompson  _et al._ (2022) A. P. Thompson, H. M. Aktulga, R. Berger, D. S. Bolintineanu, W. M. Brown, P. S. Crozier, P. J. in ’t Veld, A. Kohlmeyer, S. G. Moore, T. D. Nguyen, R. Shan, M. J. Stevens, J. Tranchida, C. Trott, and S. J. Plimpton, “LAMMPS - a flexible simulation tool for particle-based materials modeling at the atomic, meso, and continuum scales,” Comput. Phys. Commun. 271, 108171 (2022).
  * Frenkel and Smit (2002) D. Frenkel and B. Smit,  _Understanding Molecular Simulation_ (Academic Press, London, 2002).
  * Hoover (1985) W. G. Hoover, “Canonical dynamics: Equilibrium phase-space distributions,” Phys. Rev. A 31, 1695–1697 (1985).
  * Nosé (1984) S. Nosé, “A unified formulation of the constant temperature molecular dynamics methods,” J. Chem. Phys. 81, 511–519 (1984).
  * Andersen (1980) H. C. Andersen, “Molecular dynamics simulations at constant pressure and/or temperature,” J. Chem. Phys. 72, 2384–2393 (1980).
  * Kataoka and Yamada (2014) Y. Kataoka and Y. Yamada, “Phase diagram of a Lennard-Jones system by molecular dynamics simulations,” J. Comput. Chem. Jpn. 13, 115–123 (2014).
  * Travesset (2014) A. Travesset, “Phase diagram of power law and Lennard-Jones systems: Crystal phases,” J. Chem. Phys. 141, 164501 (2014).
  * Schultz and Kofke (2018) A. J. Schultz and D. A. Kofke, “Comprehensive high-precision high-accuracy equation of state and coexistence properties for classical Lennard-Jones crystals and low-temperature fluid phases,” J. Chem. Phys. 149, 204508 (2018).
  * Stephan, Staubach, and Hasse (2020) S. Stephan, J. Staubach, and H. Hasse, “Review and comparison of equations of state for the Lennard-Jones fluid,” Fluid Phase Equilibria 523, 112772 (2020).
  * Stephan  _et al._ (2019) S. Stephan, M. Thol, J. Vrabec, and H. Hasse, “Thermophysical properties of the Lennard-Jones fluid: Database and data assessment,” J. Chem. Inf. Model. 59, 4248–4265 (2019).
  * Baidakov and Tipeev (2012) V. G. Baidakov and A. O. Tipeev, “Crystal nucleation and the solid–liquid interfacial free energy,” J. Chem. Phys. 136, 074510 (2012).
  * Malek, Morrow, and Saika-Voivod (2015) S. M. A. Malek, G. P. Morrow, and I. Saika-Voivod, “Crystallization of Lennard-Jones nanodroplets: From near melting to deeply supercooled,” J. Chem. Phys. 142, 124506 (2015).
  * Morris and Song (2003) J. R. Morris and X. Song, “The anisotropic free energy of the Lennard-Jones crystal-melt interface,” J. Chem. Phys. 119, 3920–3925 (2003).
  * Larsen, Schmidt, and Schiøtz (2016) P. M. Larsen, S. Schmidt, and J. Schiøtz, “Robust structural identification via polyhedral template matching,” Modelling Simul. Mater. Sci. Eng. 24, 055007 (2016).
  * Chung  _et al._ (2022) H. W. Chung, R. Freitas, G. Cheon, and E. J. Reed, “Data-centric framework for crystal structure identification in atomistic simulations using machine learning,” Phys. Rev. Materials 6, 043801 (2022).
  * Çamkıran, Parsch, and Hibbard (2022) J. Çamkıran, F. Parsch, and G. D. Hibbard, “A local orientational order parameter for systems of interacting particles,” J. Chem. Phys. 156, 091101 (2022).
  * Stukowski (2010) A. Stukowski, “Visualization and analysis of atomistic simulation data with OVITO–the Open Visualization Tool,” Modelling Simul. Mater. Sci. Eng. 18, 015012 (2010).
  * Bulutoglu  _et al._ (2023) P. S. Bulutoglu, A. S. Zalte, N. K. Nere, D. Ramkrishna, and D. S. Corti, “A comprehensive modeling approach for polymorph selection in Lennard-Jones crystallization,” J. Chem. Phys. 158, 134505 (2023).
  * Desgranges and Delhommelle (2007) C. Desgranges and J. Delhommelle, “Controlling polymorphism during the crystallization of an atomic fluid,” Phys. Rev. Lett. 98, 235502 (2007).
  * Leyssale, Delhommelle, and Millot (2005) J.-M. Leyssale, J. Delhommelle, and C. Millot, “Atomistic simulation of the homogeneous nucleation and of the growth of N2 crystallites,” J. Chem. Phys. 122, 104510 (2005).
  * Blow  _et al._ (2023) K. E. Blow, G. A. Tribello, G. C. Sosso, and D. Quigley, “Interplay of multiple clusters and initial interface positioning for forward flux sampling simulations of crystal nucleation,” J. Chem. Phys. 158, 224102 (2023).


C_{2}(1-\gamma)\cdot\kappa(B^{\pi^{*}})\cdot\|\bm{V}_{r}(\pi^{*})\|_{\infty}\cdot\left(\frac{\|\Delta
B\|_{\infty}}{\|B^{\pi^{*}}\|_{\infty}}+\frac{\|\Delta\hat{\bm{r}}\|_{\infty}}{\|\hat{\bm{r}}(\pi^{*})\|_{\infty}}\right).$

It is clear to see that

$\|\bm{V}_{r}(\pi^{*})\|_{\infty}\leq\frac{1}{1-\gamma}$ (143)

and

$\|\bm{V}_{r}(\pi^{*})\|_{\infty}\leq\frac{\|\hat{\bm{r}}(\pi^{*})\|_{\infty}}{1-\gamma}.$
(144)

Following Jin and Sidford [2020], we have the following bound.

$\sigma^{\prime}=\|(B^{\pi^{*}})^{-1}\|_{\infty}\leq\frac{1}{1-\gamma}.$ (145)

Also, from the definition, we have that

$\|B^{\pi^{*}}\|_{\infty}=1-\gamma.$ (146)

Plugging the bound (138), (139), (143), and (144), into the inequality (142),
we have that

$\text{Regret}_{r}(\bar{\pi}^{N},N)\leq
C_{3}\cdot\frac{(|\mathcal{S}|+K)^{3}}{\alpha^{2}\cdot\xi\cdot\sigma(1-\gamma)\cdot\min\\{\sigma^{2},(1-\gamma)^{2}\cdot\delta^{2}\\}}\cdot\frac{\log^{2}(N)}{N}$
(147)

where $C_{3}$ is a constant. In a same way, for any $k\in[K]$, we obtain that

$\text{Regret}_{k}(\bar{\pi}^{N},N)\leq
C_{3}\cdot\frac{(|\mathcal{S}|+K)^{3}}{\alpha^{2}\cdot\xi\cdot\sigma(1-\gamma)\cdot\min\\{\sigma^{2},(1-\gamma)^{2}\cdot\delta^{2}\\}}\cdot\frac{\log^{2}(N)}{N}.$
(148)

To show the sample complexity bound, we denote by $\varepsilon$ a constant
such that

$\varepsilon=C_{3}\cdot\frac{(|\mathcal{S}|+K)^{3}}{\alpha^{2}\cdot\xi\cdot\sigma(1-\gamma)\cdot\min\\{\sigma^{2},(1-\gamma)^{2}\cdot\delta^{2}\\}}\cdot\frac{\log^{2}(N)}{N}.$

Therefore, we have

$N=O\left(\frac{(|\mathcal{S}|+K)^{3}}{\alpha^{2}\cdot\xi\cdot\sigma(1-\gamma)\cdot\min\\{\sigma^{2},(1-\gamma)^{2}\cdot\delta^{2}\\}}\cdot\frac{\log^{2}(1/\varepsilon)}{\varepsilon}\right).$

Note that in each of the $N$ rounds, we obtain a sample for each
$(s,a)\in\mathcal{S}\times\mathcal{A}$. Therefore, the bound on $N$ above
should be multiplied by $|\mathcal{S}|\cdot|\mathcal{A}|$ to obtain the final
sample complexity bound. Our proof is thus completed.


# Searches for BSM physics through CP violation at CDF.

, on behalf of the CDF collaboration
University and INFN Pisa
E-mail:

###### Abstract:

The CDF experiment at the Tevatron $p\bar{p}$ collider has pioneered and
established the role of hadron collisions in exploring flavor physics through
a broad program that is now at its full maturity. We report new results
sensitive to physics beyond the standard model, obtained using the whole CDF
data set; including new bounds on the mixing phase and the decay width
difference of mass-eigenstates and a measurement of the difference of CP
asymmetries in $K^{+}K^{-}$ and $\pi^{+}\pi^{-}$ decays of $D^{0}$ mesons. We
also present a new measurement of the $\Bs\rightarrow
D^{(*)+}_{s}D^{(*)-}_{s}$ branching ratio using 6.8 fb-1 of data and search
for CP violation in $D^{0}\rightarrow K^{0}_{s}\pi^{+}\pi^{-}$ decays in 6.0
fb-1 of data.

###### pacs:

13.25.Hw 11.30.Er 14.40.Nd 14.65.Fy 12.38.Qk 13.25.Ft 14.40.Lb 12.15.Ff

##  1 Introduction

Flavor physics of quarks provides one of the most promising probes for
indirect signs of new particles or interactions beyond the standard model
(SM). Strange and charm meson dynamics in particular offers rich opportunities
since its experimental exploration has not reached in extension and precision
that has been achieved in processes involving charged and neutral kaons and
bottom mesons.

The $B^{0}_{s}$ oscillations are explained in terms of second-order weak
processes involving the CKM matrix element $V_{ts}$. A broad class of generic
extensions of the SM is expected to affect the mixing amplitude, modifying the
mixing “intensity”–that is the oscillation frequency–and the phase,
$\beta_{s}=arg[-(V_{ts}V^{*}_{tb})/(V_{cs}V^{*}_{cb})]$. A non-SM enhancement
of $\beta_{s}$ would also decrease the size of the decay-width difference
$\Delta\Gamma_{s}$ between the light and heavy mass eigenstates of the
$B^{0}_{s}$ meson [1]. While the oscillation intensity has been measured
precisely [2], only loose constraints on the phase and width difference were
available until recently. The most effective determination of $\beta_{s}$ and
$\Delta\Gamma_{s}$ is achieved through the analysis of the time evolution of
“flavor–tagged” $B^{0}_{s}\kern-3.00003pt\to\kern-1.99997ptJ\\!/\\!\psi\mbox{
}\phi$ decays. The first such analysis was performed by CDF in 2008 [3]. D0
followed with a similar measurement [4] soon after CDF. In 2010 the
combination of CDF and D0 results suggested a mild deviation from the SM
expectation. However, updated measurements [4, 5, 6] showed increased
consistency with the SM, calling for additional experimental information to
clarify the picture. Here we report the latest CDF update using the final data
set of 10 fb-1. We also report a measurement of the $B^{0}_{s}\rightarrow
D^{(*)+}_{s}D^{(*)-}_{s}$ branching ratios using 6.8 fb-1 of CDF data, which
provides information on the “width–difference”.

The observation of a sizable $D^{0}$–$\bar{D}^{0}$ mixing [7, 8, 9] has raised
an increasing interest in charm dynamics, where CP violation may play an
important role. Within the SM, “CP–violating” effects are predicted to be
small since the charm transitions are described, to an excellent
approximation, by physics of the first two generations of quarks. The size of
CP violation expected from the Cabibbo-Kobayashi-Maskawa hierarchy is
$\mathcal{O}(10^{-3})$ or less. However both the $D^{0}$–$\bar{D}^{0}$ mixing
amplitude and the SM-suppressed penguin amplitude can be greatly enhanced by
new dynamics, which can also increase the size of the CP violation.

Among the most sensitive probes of physics beyond the SM are the tree-
dominated decays $D^{0}\to\pi^{+}\pi^{-}$ and $D^{0}\to K^{+}K^{-}$. Any
measured direct asymmetry significantly larger than 1% would be indication of
NP. Last year, using 5.9 fb-1 of data, CDF produced the world’s most precise
measurements of the CP asymmetries $A_{\rm{CP}}(KK)=(-0.24\pm 0.22\pm 0.09)\%$
and $A_{\rm{CP}}(\pi\pi)=(0.22\pm 0.24\pm 0.11)\%$ [10]. In spite of the
hadronic uncertainties, there is some consensus that direct CP asymmetries of
$D^{0}\to K^{+}K^{-}$ and of $D^{0}\to\pi^{+}\pi^{-}$ should be of opposite
sign. Therefore, a measurement of the difference between asymmetries of those
decays is maximally sensitive to detect direct CP violation. Indeed, the LHCb
collaboration reported recently the first evidence of CP violation in charm
measuring $\Delta A_{\rm{CP}}=A_{\rm{CP}}(KK)-A_{\rm{CP}}(\pi\pi)=(-0.82\pm
0.21\pm 0.11)\%$ [11]. An independent measurement is crucial to establish the
effect, and the 10 fb-1 sample of hadronic $D$ decays collected by CDF is the
only one currently available to attain sufficient precision.

Another suitable channel to search for CP violation in charm is the
“three–body” decay $D^{0}\rightarrow K^{0}_{s}\pi^{+}\pi^{-}$ and its resonant
substructure. In sec.6 we report the first such measurement in hadron
collisions using 6 fb-1 of data.

##  2 The CDF II detector

The collider detector at fermilab is a multipurpose experiment designed to
study $\sqrt{s}=$1.96TeV $p\overline{p}$ collisions produced by the Tevatron
collider in 2001–2011. Among the components and capabilities of the CDF II
detector [13], the tracking system is the one most relevant to these analyses.
It lies within a uniform, axial magnetic field of 1.4 T strength. The inner
tracking volume up to a radius of 28 cm is composed of 6-7 layers of
“double–sided” silicon “micro–strip” detectors [14]. An additional layer of
“single–sided” silicon is mounted directly on the “beam-pipe” at a radius of
1.5 cm, allowing excellent resolution on the impact parameter $d_{0}$, defined
as the distance of closest approach of a reconstructe “charged–particle” track
to the interaction point in the plane transverse to the beam line. The silicon
detector provides a vertex resolution of approximately 15 $\mu m$ in the
transverse and 70 $\mu m$ in the longitudinal direction. The remainder of the
tracking volume from a radius of 40 to 137 cm is occupied by an “open–cell”
drift chamber (COT) [15], providing a transverse momentum resolution of
$\sigma_{p_{T}}/p_{T}\approx 0.07\%p_{T}$ ($p_{T}$ in GeV/$c$). Hadron
identification, which is crucial for distinguishing slow kaons and protons
from pions and muons, is achieved by a likelihood combination of information
from a time-of-flight system [16] and ionization energy loss in the COT. This
offers about 1.5$\sigma$ separation between kaons and pions. A “three–level”
trigger system is used for the online event selection. At level 1 the most
important device for the described analyses is the extremely fast tracker
(XFT) [17]. It identifies charged particles using information from the COT and
measures their transverse momenta and azimuthal angles around the beam
direction. The basic requirement at level 1 is two charged particles with
transverse momentum greater than 1.5 GeV/$c$. At level 2 and 3, different
requirements are imposed for the different analyses: two “oppositely–charged”
particles reconstructed in the “drift–chamber” matched to “muon–chamber” track
segments with a dimuon mass consistent with the $J/\psi$ mass (“low–$p_{T}$”
dimuon trigger) for the analysis in sec. 3; for the remaining analyses, the
trigger uses XFT tracks combined with tracks reconstructed in the silicon
tracker, thus allowing the precise measurement of impact parameters of tracks
[18]. Impact parameters are required to be between 0.1 and 1 mm and to be
consistent with coming from a common vertex displaced from the interaction
point by at least 100 $\mu$m in the plane transverse to the beam line.

##  3 Measurement of the
$B^{0}_{s}\kern-3.00003pt\to\kern-1.99997ptJ\\!/\\!\psi\mbox{ }\phi$ time-
evolution in the final CDF Run II data set

The $B^{0}_{s}\kern-3.00003pt\to\kern-1.99997ptJ\\!/\\!\psi\mbox{ }\phi$
decays are fully reconstructed using four tracks originating from a common
displaced vertex, two matched to muon pairs consistent with a $J\\!/\\!\psi$
decay ($3.04<m_{\mu\mu}<3.14$ $\mathrm{Ge\kern-1.00006ptV\\!/}c^{2}$), and two
consistent with a $\phi\kern-3.00003pt\to\kern-1.99997ptK^{+}K^{-}$ decay
($1.009<m_{KK}<1.028$ $\mathrm{Ge\kern-1.00006ptV\\!/}c^{2}$). The dimuon mass
constraint to the known $J\\!/\\!\psi$ mass, combined with the good $p_{T}$
resolution, yield a mass resolution of the signals of about 9 MeV/$c^{2}$. The
$J\\!/\\!\psi K^{+}K^{-}$ mass distribution (fig. 1, left), shows a signal of
approximately 11 000 decays, overlapping a constant background dominated by
the prompt combinatorial component and smaller contributions from mis-
reconstructed $B$ decays.

Figure 1: Left: Distribution of $J\\!/\\!\psi K^{+}K^{-}$ mass with fit
projection overlaid. Right: Confidence regions at the 68% and 95% CL in the
($\beta_{s}$, $\Delta\Gamma_{s}$) plane.

The analysis relies on a joint fit to the time evolution of $B^{0}_{s}$ mesons
that resolves the fast oscillations by exploiting the 90 fs time resolution of
the CDF silicon detector for these final states. Because the $B^{0}_{s}$ meson
has spin zero and $J\\!/\\!\psi$ and $\phi$ have spin one, the
$B^{0}_{s}\kern-3.00003pt\to\kern-1.99997ptJ\\!/\\!\psi\mbox{ }\phi$ decay
involves three independent amplitudes, each corresponding to one possible
angular momentum state of the $J\\!/\\!\psi\phi$ system, which is also a
$C\\!P$-odd or $C\\!P$-even eigenstate. To enhance the sensitivity to
$\beta_{s}$, the time-evolution of the three decay amplitudes is fit
independently by exploiting differences in the distribution of the kaon’ and
muon’ decay angles. Sensitivity to $\beta_{s}$ can be further enhanced by
accounting for the difference in the time evolution of initially produced
$B^{0}_{s}$ and $\bar{B}^{0}_{s}$ mesons. The flavor of the meson at the time
of production is inferred by two independent classes of algorithms: the
opposite-side flavor tag (OST) and the same-side kaon tag (SSKT) [2]. OST
algorithms infer the initial flavor of the meson candidate from the decay
products of the $b$ hadron produced by the other $b$ quark in the event; SSKT
algorithms deduce the production flavor exploiting the “charge–flavor”
correlation of the neighboring kaons produced in the fragmentation process.
The OST performance has been determined with 82 000
$B^{\pm}\kern-3.00003pt\to\kern-1.99997ptJ\\!/\\!\psi(\kern-3.00003pt\to\kern-1.99997pt\mu^{+}\mu^{-})K^{\pm}$
decays fully reconstructed in the same sample as the signal. We found an
efficiency of $\varepsilon_{\rm{OST}}=(92.8\pm 0.1)\%$, an observed averaged
dilution, $D_{\rm{OST}}=1-2w_{\rm{OST}}$, equal to $(12.3\pm 0.4)\%$ and a
resulting effective tagging power of
$\varepsilon_{\rm{OST}}D^{2}_{\rm{OST}}=(1.39\pm 0.05)\%$. The SSKT algorithms
tag a smaller fraction of candidates with better precision. Its performance
has been previously determined [5] to be $\varepsilon_{\rm{SSKT}}=(52.2\pm
0.7)\%$, $D_{\rm{SSKT}}=(21.8\pm 0.3)\%$ and
$\varepsilon_{\rm{SSKT}}D^{2}_{\rm{SSKT}}=(3.2\pm 1.4)\%$. Since the SSKT
algorithm has been calibrated for early data only, we conservatively restrict
its use to the events collected in that period. Simulation shows that this
results in a modest degradation in $\beta_{s}$ resolution.

The unbinned maximum likelihood joint fit uses 9 observables from each event
to determine 32 parameters including $\beta_{s}$ and $\Delta\Gamma$, other
physics parameters ($B^{0}_{s}$ lifetime, decay amplitudes at $t=0$ and
phases, etc), and several other (“nuisance”) parameters (experimental scale
factors, etc.). If $\beta_{s}$ is fixed to its SM value, the fit shows
unbiased estimates and Gaussian uncertainties for $\Delta\Gamma_{s}$ and
$\tau_{s}$. We found $\Delta\Gamma_{s}=0.068\pm 0.026(stat)\pm 0.007(syst)$
ps-1, and mean $B^{0}_{s}$ lifetime, $\tau_{s}=1.528\pm 0.019(stat)\pm
0.009(syst)$ ps. Systematic uncertainties include mismodeling of the signal
mass model, lifetime resolution, acceptance description, and angular
distribution of the background; a $\mathcal{O}(2\%)$ contamination by
$B^{0}\kern-3.00003pt\to\kern-1.99997ptJ\\!/\\!\psi K^{*}(892)^{0}$ decays
misreconstructed as
$B^{0}_{s}\kern-3.00003pt\to\kern-1.99997ptJ\\!/\\!\psi\mbox{ }\phi$ decays;
and the silicon detector misalignment. These results are among the most
precise from a single experiment. If $\beta_{s}$ is free to float in the fit,
tests in statistical trials show that the maximum likelihood estimate is
biased for the parameters of interest, and the biases depend on the true
values of the parameters. Hence, we determine confidence regions in the
$\beta_{s}$ and $(\beta_{s},\Delta\Gamma_{s})$ spaces (fig. 1, right), by
using a profile-likelihood ratio statistic as a $\chi^{2}$ variable and
considering all other likelihood variables as nuisance parameters. Confidence
regions are corrected for “non–Gaussian” tails and systematic uncertainties to
ensure nominal coverage. By treating $\Delta\Gamma_{s}$ as a nuisance
parameter, we also obtain
$\beta_{s}\in[-\pi/2,-1.51]\cup[-0.06,0.30]\cup[1.26,\pi/2]$ rad at the 68%
CL, and $\beta_{s}\in[-\pi/2,-1.36]\cup[-0.21,0.53]\cup[1.04,\pi/2]$ rad at
the 95% CL. The fit also includes the CP-odd component that can originate by
non resonant $K^{+}K^{-}$ pair or by the $f_{0}(980)$ decays. The resulting
S-wave decay amplitude is found to be negligible. All results are consistent
with the SM expectation and with determinations of the same quantities from
other experiments [4, 6, 19].

##  4 Measurement of $B^{0}_{s}\rightarrow D^{(*)+}_{s}D^{(*)-}_{s}$ branching
ratios

A measurement of $B^{0}_{s}$ production rate times the $B^{0}_{s}\rightarrow
D^{(*)+}_{s}D^{(*)-}_{s}$ branching ratio relative to the normalization mode
$B^{0}\rightarrow D^{+}_{s}D^{-}$ is perfomed using a data sample
corresponding to an integrated luminosity of 6.8 fb-1 recorded by the
displaced track trigger [20]. $D^{+}_{s}\rightarrow K^{+}K^{-}\pi^{+}$ and
$D^{+}\rightarrow K^{-}\pi^{+}\pi^{+}$ decays are reconstructed from
combinations of three tracks with appropriate charge and mass hypothesis
assignments, fitted to a common vertex and then combined into another vertex
to form $B^{0}_{s}$ candidates. For the first time in this channel we exploit
the Dalitz structure of the intermediate states “three–body” decays is
exploited for an accurate evaluation of acceptances and efficiencies. The
relative branching fractions are determined in a simultaneous maximum
likelihood fit to the signal, $(\phi\pi^{+})(\phi\pi^{-})$ and
$(\bar{K}^{*0}K^{+})(\phi\pi^{-})$, and background,
$(\phi\pi^{+})(K^{+}\pi^{-}\pi^{-})$ and
$(\bar{K}^{*0}K^{+})(K^{+}\pi^{-}\pi^{-})$, unbinned mass distributions. The
components of the fit functions for each mass distribution are fully and
partially reconstructed signals, reflections, and background. Figure 2 shows
the projections of the fit overlaid to data. The statistical significance of
each signal exceeds 10$\sigma$.

Figure 2: Invariant mass distribution of $B^{0}_{s}\ \rightarrow
D^{+}_{s}(\phi\pi^{+})D^{-}_{s}(\phi\pi^{-})$, $B^{0}_{s}\ \rightarrow
D^{+}_{s}(\bar{K}^{*0}K^{+})D^{-}_{s}(\phi\pi^{-})$, $B^{0}\rightarrow
D^{+}_{s}(\phi\pi^{+})D^{-}(K^{+}\pi^{-}\pi^{-})$ and $B^{0}\rightarrow
D^{+}_{s}(\bar{K}^{*0}K^{+})D^{-}(K^{+}\pi^{-}\pi^{-})$ candidates with the
simultaneous fit projection overlaid. The broader structures stem from decays
where the photon or $\pi^{0}$ from the $D^{*+}_{(s)}$ decay is not
reconstructed.

Using measured values of production and relative branching fractions, the
following absolute branching fractions are derived:

$\displaystyle\mathcal{B}(B_{s}\rightarrow D^{+}_{s}D^{-}_{s})$
$\displaystyle=$ $\displaystyle(0.49\pm 0.006\pm 0.05\pm 0.08)\%,$
$\displaystyle\mathcal{B}(B_{s}\rightarrow D^{*\pm}_{s}D^{\mp}_{s})$
$\displaystyle=$ $\displaystyle(1.12\pm 0.12\pm 0.09\pm 0.19)\%,$
$\displaystyle\mathcal{B}(B_{s}\rightarrow D^{*+}_{s}D^{*-}_{s})$
$\displaystyle=$ $\displaystyle(1.79\pm 0.19\pm 0.27\pm 0.29)\%,$
$\displaystyle\mathcal{B}(B_{s}\rightarrow D^{(*)+}_{s}D^{(*-)}_{s})$
$\displaystyle=$ $\displaystyle(3.38\pm 0.25\pm 0.30\pm 0.56)\%,$

in which statistical, systematic and normalization uncertainties are reported.
These results are the most precise to date from a single experiment and can
provide information related to the “decay–width” difference
$\Delta\Gamma_{s}$.

##  5 Measurement of CP violation in charm decays in the final CDF Run II data
set

CDF previously measured CP violation in $D^{0}\to K^{+}K^{-}$ and
$D^{0}\to\pi^{+}\pi^{-}$ decays. In a larger data sample, we have measured the
difference between those asymmetries with greater precision. The analysis
follows closely the measurement of individual asymmetries [21]. The flavor of
the $D^{0}$ meson is tagged from the charge of the soft pion in the strong
$D^{\star+}\to D^{0}\pi^{+}$ decay. Since $D^{\star+}$ and $D^{\star-}$ mesons
are produced in equal number in $p\bar{p}$ interactions, any asymmetry between
the number of $D^{0}$ and $\bar{D}^{0}$ decays is due to either CP violation
or instrumental effects. The latter can be induced only by the difference in
reconstruction efficiency between positive and negative soft pions. Provided
that the relevant kinematic distributions are equalized in the two decay
channels, the instrumental asymmetry cancels to an excellent level of accuracy
in the difference between the observed asymmetries between signal yields. Such
cancellation allows an increase in sensitivity on $\Delta A_{\rm{CP}}$ by
loosening some selection criteria with respect to the measurement of
individual asymmetries thus doubling the signal yields. The offline selection
follows the standard trigger selection with some basic additional requirements
on track and vertex quality. The numbers of $D^{0}$ and $\bar{D}^{0}$ decays
are determined with a simultaneous fit to the $D^{0}\pi$-mass distribution of
positive and negative $D^{\star}$ decays. About 1.21 times $10^{6}$ $D^{0}\to
K^{+}K^{-}$ decays, fig. 3 bottom, and 550 times $10^{3}$
$D^{0}\to\pi^{+}\pi^{-}$ decays, fig. 3 top, are reconstructed, yielding the
following observed asymmetries between signal yields,
$A_{\rm{raw}}(KK)=(-2.33\pm 0.14)\%$ and $A_{\rm{raw}}(\pi\pi)=(-1.71\pm
0.15)\%$. Residual systematic uncertainties total 0.10% and are driven by
differences between $D^{*}$ mass distributions associated with charm and
anticharm decays. The final result is $\Delta A_{\rm{CP}}=(-0.62\pm 0.21\pm
0.10)\%$, which is 2.7$\sigma$ different from zero [22]. This provides strong
indication of CP violation in CDF charm data, supporting the LHCb earlier
evidence with the same resolution. The combination of CDF, LHCb, and
$B$-factory measurements deviates by approximately 3.8$\sigma$ from the no CP
violation point.

Figure 3: Projections of the combined fit on data for tagged $D^{0}\rightarrow
K^{+}K^{-}$ (top) and $D^{0}\rightarrow\pi^{+}\pi^{-}$ (bottom) decays. Charm
decays on the left and anticharm on the right.

##  6 Search for CP violation in $D^{0}\rightarrow K^{0}_{s}\pi^{+}\pi^{-}$
decays

We exploit a large sample of $D^{*+}$ decays with the $D^{0}$ production
flavor determined by the charge of the soft pion in $D^{*+}\rightarrow
D^{0}\pi^{+}$ decays. The offline reconstruction of candidates starts with
refitting tracks using the pion mass hypothesis. Two oppositely charged tracks
are combined to form a $K^{0}_{S}$ candidate. To construct $D^{0}$ candidates,
each $K^{0}_{S}$ candidate is then combined with all possible
“oppositely–charged track–pairs” from the remaining tracks in the event.
Finally, the $D^{*+}$ candidates are obtained by combining each $D^{0}$
candidate with one of the still remaining tracks in the event. The tracks
forming the $K^{0}_{S}$ , $D^{0}$ , and $D^{*+}$ candidates are subjected to
kinematic fits that constrain them to originate from common vertices. Standard
quality requirements on tracks and vertices are used to ensure “well–measured”
masses and “decay–positions”. An appropriately trained neural network
classifer contributes the final discrimination between signal and background.

The resonant substructure of a three-body decay is described using the Dalitz
plot method [23]. The Dalitz plot of the considered decay $D^{0}\rightarrow
K^{0}_{s}\pi^{+}\pi^{-}$, composed of all selected candidates, is shown in
fig. 4. Three types of intermediate resonances contribute: Cabibbo allowed,
doubly Cabibbo suppressed, and CP eigenstates. The dominant decay mode is the
Cabibbo allowed $D^{0}\rightarrow K^{*}(892)^{-}\pi^{+}$ which amounts to
about 60% of the total branching fraction. The second largest contribution is
from the intermediate CP eigenstate $K^{0}_{S}\rho(770)$, which is color
suppressed compared to $K^{∗}(892)^{-}\pi^{+}$.

Figure 4: Dalitz plot of the decay $D^{0}\rightarrow K^{0}_{s}\pi^{+}\pi^{-}$
, where the squared invariant masses of the “two–body” combinations
$K^{0}_{S}\pi^{-}$ and $\pi^{+}\pi^{-}$ are chosen as kinematic quantities.
The red line indicates the kinematic boundaries.

A binned maximum likelihood fit to the “two–dimensional” Dalitz plot
distribution with bin widths of 0.025 GeV2/c4 in both dimensions is performed
to determine the contributions of the different intermediate resonances. The
isobar model is used to describe the form of the matrix element implemented in
the likelihood. The results of the simultaneously fit to $D^{0}$ and
$\bar{D}^{0}$ Dalitz plots for the CP violating amplitudes and phases are
reported in Ref. [24]. The overall integrated CP asymmetry is found to be
$A_{CP}=(-0.05\pm 0.57\pm 0.54)\%$ with a large improvement over previous
results [25]. All the CP violating quantities are found to be consistent with
zero. Following the so-called Miranda procedure [26], we also performed a
“model–independent” search for CP violation in the Dalitz plot distribution of
the decay $D^{0}\rightarrow K^{0}_{s}\pi^{+}\pi^{-}$ by comparing the binned,
normalized Dalitz plots for $D^{0}$ and $\bar{D}^{0}$. No assumptions about
the resonant substructure of the decay are used. This approach confirms that
“no–CP–violation“ between the $D^{0}$ and $\bar{D}^{0}$ Dalitz plots is
present.

##  7 Summary

The results of a few recent updates on flagship flavor measurements that use
the complete CDF II data set or a large fraction of it have been reported.
Improved bounds on the $B^{0}_{s}$ mixing phase and decay width difference of
$B^{0}_{s}$ mass-eigenstates are found to be consistent with the standard
model. We also presented the most precise measurement of $B^{0}_{s}\rightarrow
D^{(*)+}_{s}D^{(*)-}_{s}$ branching ratio. A measurement of the difference of
CP asymmetries in $K^{+}K^{-}$ and $\pi^{+}\pi^{-}$ decays of $D^{0}$ meson
shows significant discrepancy from zero confirming and supporting similar
results from LHCb. Its interpretation as due to physics within or beyond the
SM is still under debate. We finally presented a null search for CP violation
in $D^{0}\rightarrow K^{0}_{s}\pi^{+}\pi^{-}$ decays in 6.0 fb-1 of data.

These measurements are a sampling from an intensive flavor physics program
developed at the Tevatron in the last decade that is now at its full maturity.
The chief achievement has been to extend and deepen the study of the
$B^{0}_{s}$ mesons dynamics, which was only marginally explored before.
Stringent constraints on the presence of NP in $B^{0}_{s}$ mixing has been
imposed by measuring with a good precision both its intensity and its
associated phase. In addition, favorable production cross sections together
with the CP-symmetric nature of the $p\overline{p}$ collisions allowed
exploration of CPV in the charm sector with “world–leading” precision. Another
important heritage has been to show that flavor physics is in fact possible
also at hadron colliders and provides information competitive and
complementary to the information from dedicated experiments. Many trigger
reconstruction, and data analysis techniques developed at CDF are now used as
standard by the next generation of flavor experiments at hadron colliders.

## References

  * [1] S. Faller, R. Fleischer, and T. Mannel, Phys. Rev. D 79, 014005 (2009).
  * [2] A. Abulencia and others, _CDF Collaboration_ , Phys. Rev. Lett. 97, 242003 (2006).
  * [3] T. Aaltonen and others, _CDF Collaboration_ , Phys. Rev. Lett. 100, 161802 (2008); V.M. Abazov and others, _D0 Collaboration_ , Phys. Rev. Lett. 101, 241801 (2008).
  * [4] V.M. Abazov and others, _D0 Collaboration_ , Phys. Rev. D 85, 032006 (2012), Phys. Rev. D 82, 082001 (2010), Phys. Rev. D 84, 052007 (2011).
  * [5] T. Aaltonen and others, _CDF Collaboration_ , Phys. Rev. D 85, 072002 (2012).
  * [6] R. Aaij and others, _LHCb Collaboration_ , Phys. Rev. Lett. 108, 101803 (2012).
  * [7] B. Auber and others, _BABAR Collaboration_ , Phys. Rev. Lett. 98, 211802 (2007).
  * [8] M. Staric and others, _Belle Collaboration_ , Phys. Rev. Lett. 98, 211803 (2007).
  * [9] T. Aaltonen and others, _CDF Collaboration_ , Phys. Rev. Lett. 100, 121802 (2008).
  * [10] T. Aaltonen and others, _CDF Collaboration_ , Phys. Rev. D 85, 012009 (2012).
  * [11] R. Aaij and others, _LHCb Collaboration_ , Phys. Rev. Lett. 108, 111602 (2012).
  * [12] I.I Bigi and H. Yamamoto, Phys. Lett. B 349, 363 (1995).
  * [13] D. Acosta and others, _CDF Collaboration_ , Phys. Rev. D 71, 032001 (2005).
  * [14] C.S. Hill, _CDF Collaboration_ , Nucl. Instrum. Methods A 530, 1 (2004).
  * [15] A.A. Affolder and others, _CDF Collaboration_ , Nucl. Instrum. Methods A 526, 249 (2004).
  * [16] D. Acosta and others, _CDF Collaboration_ , Nucl. Instrum. Methods A 518, 605 (2004).
  * [17] E.J. Thomson and others, IEEE Trans. Nucl. Sci. 49, 1063 (2002).
  * [18] G.Punzi and L. Ristori, Annual Review of Nuclear and Particle Science 60, 595-614 (2010).
  * [19] R. Aaij and others, _LHCb Collaboration_ , [arXiv:1202.4717 (2012)].
  * [20] T. Aaltonen and others, _CDF Collaboration_ , Phys. Rev. Lett. 108, 201801 (2012).
  * [21] T. Aaltonen and others, _CDF Collaboration_ , Phys. Rev. D 85, 012009 (2012).
  * [22] T. Aaltonen and others, _CDF Collaboration_ , [arXiv:1207.2158].
  * [23] R.H Dalitz, _Phil. Mag._ 44, 1068 (1953).
  * [24] T. Aaltonen and others, _CDF Collaboration_ , [arXiv:1207.0825].
  * [25] D.M. Asner and others, _CLEO Collaboration_ , Phys. Rev. Lett. 70, 091101 (2004).
  * [26] I. Bediaga and others, Phys. Rev. D 80, 096006 (2009).